rm(list=ls())
require(dplyr)
require(CBPS)
require(AER)
require(ggplot2)
require(multiwayvcov)
require(stargazer)
require(reshape2)
require(car)
require(broom)
require(dotwhisker)
require(gridExtra)
require(gtable)
require(grid)

load("rep_data.RData")

######################
##MANUSCRIPT FIGURES##
######################

##FIGURE ONE: CREATED MECHANICALLY##

##FIGURE TWO##
bins <- analysis
bins$distance_bins <- cut(bins$distance, breaks = 19.9:61, right=T)
binner <- NULL
binner$mean_treat <- summarise(group_by(bins, distance_bins), mean(treat))
binner <- as.data.frame(binner)
binner$distance <- c(20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 45, 46, 47, 48, 49, 50, 51, 52, 53, 
                     54, 55, 56, 57, 58, 59, 60, NA) 


RD_g <- ggplot(analysis, aes(x = distance, y = treat, color = as.factor(forcing)))+  geom_jitter(height = .01, alpha = .2) +
  xlab("Distance From Pole")+ylab("Share With Legal Connection")+
  geom_vline(xintercept=40, linetype="longdash")+
  geom_point(data = binner, aes(x = distance, y = mean_treat.mean.treat.), shape=15, color = "red4")+
  scale_colour_discrete(name="Experimental\nCondition",breaks=c("0", "1"), labels=c("Control", "Forcing"))

RD_g


##FIGURE 3##
#scaled data set
analysis_scale <- dplyr::select(analysis, -gender, -birthplace, -age, -forcing, -treat, -religion, -caste, -Pole, -education, -village)
analysis_scale <- as.data.frame(lapply(analysis_scale, scale))
analysis_scale$treat <- analysis$treat
analysis_scale$gender <- analysis$gender
analysis_scale$birthplace <- analysis$birthplace
analysis_scale$age <- analysis$age
analysis_scale$forcing <- analysis$forcing
analysis_scale$religion <- analysis$religion
analysis_scale$caste <- analysis$caste
analysis_scale$Pole <- analysis$Pole

#Expenditures#
e_fit1s <- ivreg(total_expenditure ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis_scale)
e_se1s <- coeftest(e_fit1s, cluster.vcov(e_fit1s, as.factor(analysis_scale$Pole)))

#Food
e_fit2s <- ivreg(food_expenditure ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis_scale)
e_se2s <- coeftest(e_fit2s, cluster.vcov(e_fit2s, as.factor(analysis_scale$Pole)))

#Education
e_fit3s <- ivreg(education_expenditure ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis_scale)
e_se3s <- coeftest(e_fit3s, cluster.vcov(e_fit3s, as.factor(analysis_scale$Pole)))

#kerosene
e_fit4s <- ivreg(kerosene_expenditure ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis_scale)
e_se4s <- coeftest(e_fit4s, cluster.vcov(e_fit4s, as.factor(analysis_scale$Pole)))

##create objects for dwplot
ex1 <- tidy(e_fit1s)
ex1$std.error <- e_se1s[,2]
ex2 <- tidy(e_fit2s)
ex2$std.error <- e_se2s[,2]
ex3 <- tidy(e_fit3s)
ex3$std.error <- e_se3s[,2]
ex4 <- tidy(e_fit4s)
ex4$std.error <- e_se4s[,2]

#the tidy dataframe
coefex <- rbind(ex1[2,], ex2[2,], ex3[2,], ex4[2,])
coefex$term <- c("Total Expenditure", "Food Expenditure", "Education Expenditure", "Kerosene Expenditure")

lim1 <- dwplot(coefex)+geom_vline(xintercept = 0, colour = "grey50", linetype = 2)+xlab("Standard Score") + ggtitle("(A) Expenditure")+xlim(-2, 2)



#Household Activity#
#adult activity
ha_fit1s <- ivreg(adult_activity ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis_scale)
ha_se1s <- coeftest(ha_fit1s, cluster.vcov(ha_fit1s, as.factor(analysis_scale$Pole)))

#child activity
ha_fit2s <- ivreg(child_activity ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis_scale)
ha_se2s <- coeftest(ha_fit2s, cluster.vcov(ha_fit2s, as.factor(analysis_scale$Pole)))


has1 <- tidy(ha_fit1s)
has1$std.error <- ha_se1s[,2]

has2 <- tidy(ha_fit2s)
has2$std.error <- ha_se2s[,2]

coefhas <- rbind(has1[2,], has2[2,])
coefhas$term <- c("Adult Activity", "Child Activity")

lim2 <- dwplot(coefhas)+geom_vline(xintercept = 0, colour = "grey50", linetype = 2)+xlab("Standard Score")  + ggtitle("(B) Activity")+xlim(-2, 2)

#Appliances#
a_fit1s <- ivreg(appliances ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis_scale)
a_se1s <- coeftest(a_fit1s, cluster.vcov(a_fit1s, as.factor(analysis_scale$Pole)))

#appliance use
a_fit2s <- ivreg(appliance_use ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis_scale)
a_se2s <- coeftest(a_fit2s, cluster.vcov(a_fit2s, as.factor(analysis_scale$Pole)))

as1 <- tidy(a_fit1s)
as1$std.error <- a_se1s[,2]

as2 <- tidy(a_fit2s)
as2$std.error <- a_se2s[,2]

coefas <- rbind(as1[2,], as2[2,])
coefas$term <- c("Number of Appliances", "Appliance Use")

lim3 <- dwplot(coefas)+geom_vline(xintercept = 0, colour = "grey50", linetype = 2)+xlab("Standard Score") + ggtitle("(C) Appliances")+xlim(-2, 2)


lim1 <- ggplotGrob(lim1)
lim2 <- ggplotGrob(lim2)
lim3 <- ggplotGrob(lim3)
glim <- rbind(lim1, lim2, lim3, size="first")
glim$widths <-  unit.pmax(lim1$widths, lim2$widths, lim3$widths)


grid.arrange(glim, ncol=1)


######################
##MANUSCRIPT TABLES##
######################

##TABLE ONE: CREATED MECHANICALLY##


##TABLE TWO: CREATED MECHANICALLY##


##TABLE THREE## Line added in LaTeX
b1 <- tidy(t.test(gender ~ forcing, data=analysis))
b2 <- tidy(t.test(birthplace ~ forcing, data=analysis))
b3 <- tidy(t.test(age ~ forcing, data=analysis))
b4 <- tidy(t.test(religion ~ forcing, data=analysis))
b5 <- tidy(t.test(caste ~ forcing, data=analysis))
b6 <- tidy(t.test(education ~ forcing, data = analysis))
b7 <- tidy(t.test(homeBuilt ~ forcing, data = analysis))

balance_t <- rbind(b1, b2, b3, b4, b5, b6, b7)
balance_t <- balance_t[, c("estimate1", "estimate2", "p.value")]


row.names(balance_t) <- c("Female Household Head", "Birthplace", "Age", "Hindu", "Scheduled Caste/Tribe", "Education Level", "Age of Home")


stargazer(balance_t, type="latex", float = F, summary= FALSE,
          covariate.labels = c("Covariate", "Control Mean", "Treatment Mean", "p.value"))


##TABLE FOUR##
#Total Expenditure
e_fit1 <- ivreg(total_expenditure ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
e_se1 <- coeftest(e_fit1, cluster.vcov(e_fit1, as.factor(analysis$Pole)))

#Food
e_fit2 <- ivreg(food_expenditure ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
e_se2 <- coeftest(e_fit2, cluster.vcov(e_fit2, as.factor(analysis$Pole)))

#Education
e_fit3 <- ivreg(education_expenditure ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
e_se3 <- coeftest(e_fit3, cluster.vcov(e_fit3, as.factor(analysis$Pole)))

#kerosene
e_fit4 <- ivreg(kerosene_expenditure ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
e_se4 <- coeftest(e_fit4, cluster.vcov(e_fit4, as.factor(analysis$Pole)))

#electricity
e_fit5 <- ivreg(electricity_expenditure ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
e_se5 <- coeftest(e_fit5, cluster.vcov(e_fit5, as.factor(analysis$Pole)))


##Adjust p values for multiple comparison
e_pvalues <- cbind(e_se1[,4], e_se2[,4], e_se3[,4], e_se4[,4])

#p.adjust(e_pvalues[1,], method="BH", n=length(e_pvalues[1,]))

e_pvalues_adj <- t(apply(e_pvalues, MARGIN=1, function(x) p.adjust(x, method="BH", n=length(x))))


##Create the lists for stargazer
fit.list.e <- list(e_fit1, e_fit2, e_fit3, e_fit4)
se.list.e <- list(e_se1[,2], e_se2[,2], e_se3[,2], e_se4[,2])
p.list.e <- list(e_pvalues_adj[,1], e_pvalues_adj[,2], e_pvalues_adj[,3], e_pvalues_adj[,4])
depvar.list.e <- c("Total Expenditure", "Food Expenditure", "Education Expenditure", "Kerosene Expenditure")
omit.stats <- c("rsq","ser","f")
cov.list <- c("Legal Electrification")

stargazer(title = "Legal Electrification and Expenditure",
          se = se.list.e,
          header = F,
          fit.list.e,
          p = p.list.e,
          covariate.labels = cov.list,
          dep.var.labels = depvar.list.e,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          model.names = F,
          omit = c("gender","birthplace","age","religion","caste","as.factor"),
          add.lines = list(c("Control Variables", "Yes", "Yes", "Yes", "Yes"),
                           c("N Poles", "152", "152", "152", "152")),
          notes = "Standard errors clustered at pole level. Pole fixed effects included.")



##TABLE FIVE##
#adult activity
ha_fit1 <- ivreg(adult_activity ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
ha_se1 <- coeftest(ha_fit1, cluster.vcov(ha_fit1, as.factor(analysis$Pole)))

#child activity
ha_fit2 <- ivreg(child_activity ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
ha_se2 <- coeftest(ha_fit2, cluster.vcov(ha_fit2, as.factor(analysis$Pole)))

##Adjust p values for multiple comparison
ha_pvalues <- cbind(ha_se1[,4], ha_se2[,4])

ha_pvalues_adj <- t(apply(ha_pvalues, MARGIN=1, function(x) p.adjust(x, method="BH", n=length(x))))

##Create the lists for stargazer
fit.list.ha <- list(ha_fit1, ha_fit2)
se.list.ha <- list(ha_se1[,2], ha_se2[,2])
p.list.ha <- list(ha_pvalues_adj[,1], ha_pvalues_adj[,2])
depvar.list.ha <- c("Adult Activity", "Child Activity")
omit.stats <- c("rsq","ser","f")
cov.list <- c("Legal Electrification")

##Create the table
stargazer(title = "Legal Electrification and Household Activity",
          se = se.list.ha,
          header = F,
          fit.list.ha,
          p = p.list.ha,
          covariate.labels = cov.list,
          dep.var.labels = depvar.list.ha,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          model.names = F,
          omit = c("gender","birthplace","age","religion","caste","as.factor"),
          add.lines = list(c("Control Variables", "Yes", "Yes"),
                           c("N Poles", "152", "152")),
          notes = "Standard errors clustered at pole level. Pole fixed effects included."
)


##TABLE SIX##
ad1 <- ivreg(adult_home_work ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
adse1 <- coeftest(ad1, cluster.vcov(ad1, as.factor(analysis$Pole)))

ad2 <- ivreg(adult_home_cook ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
adse2 <- coeftest(ad2, cluster.vcov(ad2, as.factor(analysis$Pole)))

ad3 <- ivreg(adult_home_charge ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
adse3 <- coeftest(ad3, cluster.vcov(ad3, as.factor(analysis$Pole)))

ad4 <- ivreg(adult_home_clean ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
adse4 <- coeftest(ad4, cluster.vcov(ad4, as.factor(analysis$Pole)))

ad5 <- ivreg(adult_home_leisure ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
adse5 <- coeftest(ad5, cluster.vcov(ad5, as.factor(analysis$Pole)))

ad6 <- ivreg(adult_home_time ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
adse6 <- coeftest(ad6, cluster.vcov(ad6, as.factor(analysis$Pole)))


fit.list.ad <- list(ad1, ad2, ad3, ad4, ad5, ad6)
se.list.ad <- list(adse1[,2], adse2[,2], adse3[,2], adse4[,2], adse5[,2], adse6[,2])
depvar.list.ad <- c("Work", "Cook", "Charge", "Clean", "Leisure", "Time Spent" )
omit.stats <- c("rsq","ser","f")
cov.list <- c("Legal Electrification")

##Create the table

stargazer(title = "Legal Electrification and Household Activity Decomposition",
          se = se.list.ad,
          header = F,
          fit.list.ad,
          
          covariate.labels = cov.list,
          dep.var.labels = depvar.list.ad,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          model.names = F,
          omit = c("gender","birthplace","age","religion","caste","as.factor"),
          add.lines = list(c("Control Variables", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("N Poles", "152", "152", "152", "152", "152", "152")),
          notes = "Standard errors clustered at pole level. Pole fixed effects included."
)


##TABLE SEVEN##
a_fit1 <- ivreg(appliances ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
a_se1 <- coeftest(a_fit1, cluster.vcov(a_fit1, as.factor(analysis$Pole)))

#appliance use
a_fit2 <- ivreg(appliance_use ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
a_se2 <- coeftest(a_fit2, cluster.vcov(a_fit2, as.factor(analysis$Pole)))

##Adjust p values for multiple comparison
a_pvalues <- cbind(a_se1[,4], a_se2[,4])

#p.adjust(e_pvalues[1,], method="BH", n=length(e_pvalues[1,]))

a_pvalues_adj <- t(apply(a_pvalues, MARGIN=1, function(x) p.adjust(x, method="BH", n=length(x))))

##Create the lists for stargazer
fit.list.a <- list(a_fit1, a_fit2)
se.list.a <- list(a_se1[,2], a_se2[,2])
p.list.a <- list(a_pvalues_adj[,1], a_pvalues_adj[,2])
depvar.list.a <- c("Number of Appliances", "Appliance Use")
omit.stats <- c("rsq","ser","f")
cov.list <- c("Legal Electrification")

##Create the table
stargazer(title = "Legal Electrification and Appliances",
          se = se.list.a,
          header = F,
          fit.list.a,
          p = p.list.a,
          covariate.labels = cov.list,
          dep.var.labels = depvar.list.a,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          model.names = F,
          omit = c("gender","birthplace","age","religion","caste","as.factor"),
          add.lines = list(c("Control Variables", "Yes", "Yes"),
                           c("N Poles", "152", "152")),
           
          notes = "Standard errors clustered at pole level. Pole fixed effects included."
)


############
##APPENDIX##
############

##FIGURE A1: CREATED MECHANICALLY##

##FIGURE A2: SEE power_analysis.R##


##TABLE A1##
analysis$hb_r <- round(analysis$homeBuilt, 0)
analysis$yrs_r <- round(analysis$yrsConnected, 0)
analysis$el_r <- round(analysis$elec_hours, 0)

stargazer(analysis[c("Pole", "forcing", "treat", "katiya", "gender",
                     "birthplace", "age", "religion", "caste", "education",
                     "hb_r", "yrs_r", "el_r")],
          type="latex", float = F, 
          covariate.labels = c("Poles", "Within Legal Distance", "Legal Connection", "Illegal Connection","Female Household Head", "Birthplace", "Age", "Hindu", "Scheduled Caste/Tribe", "Education Level", "Age of House", "Years Connected to Grid", "Daily Hours of Electricity"))


##TABLE A2: CREATED MECHANICALLY##

##FIGURE A3##
##Expenditures##
e_fit1s <- ivreg(total_expenditure ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis_scale)
e_se1s <- coeftest(e_fit1s, cluster.vcov(e_fit1s, as.factor(analysis_scale$Pole)))

#Food
e_fit2s <- ivreg(food_expenditure ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis_scale)
e_se2s <- coeftest(e_fit2s, cluster.vcov(e_fit2s, as.factor(analysis_scale$Pole)))

#Education
e_fit3s <- ivreg(education_expenditure ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis_scale)
e_se3s <- coeftest(e_fit3s, cluster.vcov(e_fit3s, as.factor(analysis_scale$Pole)))

#kerosene
e_fit4s <- ivreg(kerosene_expenditure ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis_scale)
e_se4s <- coeftest(e_fit4s, cluster.vcov(e_fit4s, as.factor(analysis_scale$Pole)))

##create objects for dwplot
ex1 <- tidy(e_fit1s)
ex1$std.error <- e_se1s[,2]
ex2 <- tidy(e_fit2s)
ex2$std.error <- e_se2s[,2]
ex3 <- tidy(e_fit3s)
ex3$std.error <- e_se3s[,2]
ex4 <- tidy(e_fit4s)
ex4$std.error <- e_se4s[,2]

#the tidy dataframe
coefex <- rbind(ex1[2,], ex2[2,], ex3[2,], ex4[2,])
coefex$term <- c("Total Expenditure", "Food Expenditure", "Education Expenditure", "Kerosene Expenditure")


g1 <- dwplot(coefex)+geom_vline(xintercept = 0, colour = "grey50", linetype = 2)+xlab("Standard Score") + ggtitle("(A) Expenditure")+xlim(-2, 2)

##KEROSENE##
ke_fit1s <- ivreg(kerosene_lamps ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis_scale)
ke_se1s <- coeftest(ke_fit1s, cluster.vcov(ke_fit1s, as.factor(analysis_scale$Pole)))

#num kerosene lamps
ke_fit2s <- ivreg(num_kerosene_lamps ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis_scale)
ke_se2s <- coeftest(ke_fit2s, cluster.vcov(ke_fit2s, as.factor(analysis_scale$Pole)))

#kerosene lamp hours
ke_fit3s <- ivreg(kerosene_lamp_hours ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis_scale)
ke_se3s <- coeftest(ke_fit3s, cluster.vcov(ke_fit3s, as.factor(analysis_scale$Pole)))

#Other use of Kerosene
ke_fit4s <- ivreg(kerosene_other ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis_scale)
ke_se4s <- coeftest(ke_fit4s, cluster.vcov(ke_fit4s, as.factor(analysis_scale$Pole)))

kes1 <- tidy(ke_fit1s)
kes1$std.error <- ke_se1s[,2]

kes2 <- tidy(ke_fit2s)
kes2$std.error <- ke_se2s[,2]

kes3 <- tidy(ke_fit3s)
kes3$std.error <- ke_se3s[,2]

kes4 <- tidy(ke_fit4s)
kes4$std.error <- ke_se4s[,2]

coefke <- rbind(kes1[2,], kes2[2,], kes3[2,], kes4[2,])
coefke$term <- c("Kerosene Lamp Use", "Number of Kerosene Lamps", "Kerosene Lamp Hours", "Other Use of Kerosene")


g2 <- dwplot(coefke)+geom_vline(xintercept = 0, colour = "grey50", linetype = 2)+xlab("Standard Score") + ggtitle("(B) Kerosene")+xlim(-2, 2)

##Artificial Lighting##
#Lighting hours
l_fit1s <- ivreg(lighting_hours ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis_scale)
l_se1s <- coeftest(l_fit1s, cluster.vcov(l_fit1s, as.factor(analysis_scale$Pole)))

#child lighting hours
l_fit2s <- ivreg(child_lighting ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis_scale)
l_se2s <- coeftest(l_fit2s, cluster.vcov(l_fit2s, as.factor(analysis_scale$Pole)))

#adult lighting hours
l_fit3s <- ivreg(adult_lighting ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis_scale)
l_se3s <- coeftest(l_fit3s, cluster.vcov(l_fit3s, as.factor(analysis_scale$Pole)))



ls1 <- tidy(l_fit1s)
ls1$std.error <- l_se1s[,2]

ls2 <- tidy(l_fit2s)
ls2$std.error <- l_se2s[,2]

ls3 <- tidy(l_fit3s)
ls3$std.error <- l_se3s[,2]


coefls <- rbind(ls1[2,], ls2[2,], ls3[2,])
coefls$term <- c("Household Light", "Child Light", "Adult Light")



g3 <- dwplot(coefls)+geom_vline(xintercept = 0, colour = "grey50", linetype = 2)+xlab("Standard Score") + ggtitle("(C) Lighting")+xlim(-2, 2)


##Appliances##
a_fit1s <- ivreg(appliances ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis_scale)
a_se1s <- coeftest(a_fit1s, cluster.vcov(a_fit1s, as.factor(analysis_scale$Pole)))

#appliance use
a_fit2s <- ivreg(appliance_use ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis_scale)
a_se2s <- coeftest(a_fit2s, cluster.vcov(a_fit2s, as.factor(analysis_scale$Pole)))

as1 <- tidy(a_fit1s)
as1$std.error <- a_se1s[,2]

as2 <- tidy(a_fit2s)
as2$std.error <- a_se2s[,2]

coefas <- rbind(as1[2,], as2[2,])
coefas$term <- c("Number of Appliances", "Appliance Use")


g4 <- dwplot(coefas)+geom_vline(xintercept = 0, colour = "grey50", linetype = 2)+xlab("Standard Score") + ggtitle("(D) Appliances")+xlim(-2, 2)


##Household Activity##
#adult activity
ha_fit1s <- ivreg(adult_activity ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis_scale)
ha_se1s <- coeftest(ha_fit1s, cluster.vcov(ha_fit1s, as.factor(analysis_scale$Pole)))

#child activity
ha_fit2s <- ivreg(child_activity ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis_scale)
ha_se2s <- coeftest(ha_fit2s, cluster.vcov(ha_fit2s, as.factor(analysis_scale$Pole)))


has1 <- tidy(ha_fit1s)
has1$std.error <- ha_se1s[,2]

has2 <- tidy(ha_fit2s)
has2$std.error <- ha_se2s[,2]

coefhas <- rbind(has1[2,], has2[2,])
coefhas$term <- c("Adult Activity", "Child Activity")

g5 <- dwplot(coefhas)+geom_vline(xintercept = 0, colour = "grey50", linetype = 2)+xlab("Standard Score")  + ggtitle("(E) Activity")+xlim(-2, 2)




##Satisfaction and Attitudes##
#reliability
sa_fit1s <- ivreg(satisfaction_reliability ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis_scale)
sa_se1s <- coeftest(sa_fit1s, cluster.vcov(sa_fit1s, as.factor(analysis_scale$Pole)))

#cost
sa_fit2s <- ivreg(satisfaction_cost ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis_scale)
sa_se2s <- coeftest(sa_fit2s, cluster.vcov(sa_fit2s, as.factor(analysis_scale$Pole)))

#safety
sa_fit3s <- ivreg(satisfaction_safety ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis_scale)
sa_se3s <- coeftest(sa_fit3s, cluster.vcov(sa_fit3s, as.factor(analysis_scale$Pole)))

#brightness
sa_fit4s <- ivreg(satisfaction_brightness ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis_scale)
sa_se4s <- coeftest(sa_fit4s, cluster.vcov(sa_fit4s, as.factor(analysis_scale$Pole)))


#overall satisfaction
sa_fit5s <- ivreg(satisfaction ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis_scale)
sa_se5s <- coeftest(sa_fit5s, cluster.vcov(sa_fit5s, as.factor(analysis_scale$Pole)))

#change in satisfaction
sa_fit6s <- ivreg(satisfaction_chng ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis_scale)
sa_se6s <- coeftest(sa_fit6s, cluster.vcov(sa_fit6s, as.factor(analysis_scale$Pole)))

#electricity value
sa_fit7s <- ivreg(elec_value ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis_scale)
sa_se7s <- coeftest(sa_fit7s, cluster.vcov(sa_fit7s, as.factor(analysis_scale$Pole)))

#increase in business income
sa_fit9s <- ivreg(income_increase ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis_scale)
sa_se9s <- coeftest(sa_fit9s, cluster.vcov(sa_fit9s, as.factor(analysis_scale$Pole)))

#interest in business
sa_fit10s <- ivreg(business_interest ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis_scale)
sa_se10s <- coeftest(sa_fit10s, cluster.vcov(sa_fit10s, as.factor(analysis_scale$Pole)))

#aspirations/modernity
sa_fit11s <- ivreg(aspirations ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis_scale)
sa_se11s <- coeftest(sa_fit11s, cluster.vcov(sa_fit11s, as.factor(analysis_scale$Pole)))



sas1 <- tidy(sa_fit1s)
sas1$std.error <- sa_se1s[,2]

sas2 <- tidy(sa_fit2s)
sas2$std.error <- sa_se2s[,2]

sas3 <- tidy(sa_fit3s)
sas3$std.error <- sa_se3s[,2]

sas4 <- tidy(sa_fit4s)
sas4$std.error <- sa_se4s[,2]

sas5 <- tidy(sa_fit5s)
sas5$std.error <- sa_se5s[,2]

sas6 <- tidy(sa_fit6s)
sas6$std.error <- sa_se6s[,2]

sas7 <- tidy(sa_fit7s)
sas7$std.error <- sa_se7s[,2]


sas9 <- tidy(sa_fit9s)
sas9$std.error <- sa_se9s[,2]

sas10 <- tidy(sa_fit10s)
sas10$std.error <- sa_se10s[,2]

sas11 <- tidy(sa_fit11s)
sas11$std.error <- sa_se11s[,2]

coefsas <- rbind(sas1[2,], sas2[2,], sas3[2,], sas4[2,], sas5[2,], sas6[2,], sas7[2,], sas9[2,], sas10[2,], sas11[2,])
coefsas$term <- c("Reliability", "Cost", "Safety", "Brightness", "Overall Satisfaction", "Change in Satisfaction", "Electricity Value", "Increase in Income", "Business Aspiration", "General Aspirations")



g6 <- dwplot(coefsas)+geom_vline(xintercept = 0, colour = "grey50", linetype = 2)+xlab("Standard Score") + ggtitle("(F) Satisfaction")+xlim(-2, 2)


##Knowledge##

k_fit1s <- ivreg(knowledge ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis_scale)
k_se1s <- coeftest(k_fit1s, cluster.vcov(k_fit1s, as.factor(analysis_scale$Pole)))

ks1 <- tidy(k_fit1s)
ks1$std.error <- k_se1s[,2]

coefks <- ks1[2,]
coefks$term <- "Knowledge"


g7 <- dwplot(coefks)+geom_vline(xintercept = 0, colour = "grey50", linetype = 2)+xlab("Standard Score") + ggtitle("(G) Knowledge")+xlim(-2, 2)


# Combine plots into a single plot


p1 <- ggplotGrob(g1)
p2 <- ggplotGrob(g2)
p3 <- ggplotGrob(g3)
p4 <- ggplotGrob(g4)
p5 <- ggplotGrob(g5)
p6 <- ggplotGrob(g6)
p7 <- ggplotGrob(g7)

ga <- rbind(p1, p2, p3, p4, size="first")
ga$widths <- unit.pmax(p1$widths, p2$widths, p3$widths, p4$widths)

gb <- rbind(p5, p6, p7, size="first")
gb$widths <-  unit.pmax(p5$widths, p6$widths, p7$widths)
gb$heights <- ga$heights

grid.arrange(ga, gb, ncol=2)


##FIGURE A4##
df <- analysis[complete.cases(analysis[, 43]),]
df$fuel_type[df$fuel_type=="3"] <- "5" #collapse minigrid into other- so few observations
df2 <- df %>%
  count(fuel_type, forcing) %>%
  group_by(fuel_type) %>% 
  mutate(prop1 = n/342)
df2$prop2 <- df2$n/344

df2 <- mutate(df2, prop = ifelse(forcing==0, prop1, prop2))


fuel_labels <- c("Grid", "Kerosene", "Solar", "Other")

fuel_hist <- ggplot(data = df2, aes(as.factor(fuel_type), prop, fill = as.factor(forcing))) + 
  geom_bar(stat = "identity", position = "dodge")+xlab("Primary Fuel Source for Lighting")+ylab("Share of Households")+
  scale_fill_manual( name = "Experimental\nCondition", breaks = c("0", "1"), labels = c("Control", "Forcing"),
                     values = scales::hue_pal()(2) ) + 
  scale_x_discrete(labels= fuel_labels)



fuel_hist


##TABLE A3##

##Total expenditure##
ex_zsc_reg <- ivreg(ex_zsc ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
ex_zsc_reg_se <- coeftest(ex_zsc_reg, cluster.vcov(ex_zsc_reg, as.factor(analysis$Pole)))

ex_icw_reg <- ivreg(ex_icw ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
ex_icw_reg_se <- coeftest(ex_icw_reg, cluster.vcov(ex_icw_reg, as.factor(analysis$Pole)))

##Kerosene##
k_zsc_reg <- ivreg(k_zsc ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
k_zsc_reg_se <- coeftest(k_zsc_reg, cluster.vcov(k_zsc_reg, as.factor(analysis$Pole)))

k_icw_reg <- ivreg(k_icw ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
k_icw_reg_se <- coeftest(k_icw_reg, cluster.vcov(k_icw_reg, as.factor(analysis$Pole)))

##Artificial Lighting##
l_zsc_reg <- ivreg(l_zsc ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
l_zsc_reg_se <- coeftest(l_zsc_reg, cluster.vcov(l_zsc_reg, as.factor(analysis$Pole)))

l_icw_reg <- ivreg(l_icw ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
l_icw_reg_se <- coeftest(l_icw_reg, cluster.vcov(l_icw_reg, as.factor(analysis$Pole)))

##Appliances##
a_zsc_reg <- ivreg(a_zsc ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
a_zsc_reg_se <- coeftest(a_zsc_reg, cluster.vcov(a_zsc_reg, as.factor(analysis$Pole)))

a_icw_reg <- ivreg(a_icw ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
a_icw_reg_se <- coeftest(a_icw_reg, cluster.vcov(a_icw_reg, as.factor(analysis$Pole)))


##Household Activity##
ha_zsc_reg <- ivreg(ha_zsc ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
ha_zsc_reg_se <- coeftest(ha_zsc_reg, cluster.vcov(ha_zsc_reg, as.factor(analysis$Pole)))

ha_icw_reg <- ivreg(ha_icw ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
ha_icw_reg_se <- coeftest(ha_icw_reg, cluster.vcov(ha_icw_reg, as.factor(analysis$Pole)))


##Satisfaction and Attitudes##
sa_zsc_reg <- ivreg(sa_zsc ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
sa_zsc_reg_se <- coeftest(sa_zsc_reg, cluster.vcov(ha_zsc_reg, as.factor(analysis$Pole)))

sa_icw_reg <- ivreg(sa_icw ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
sa_icw_reg_se <- coeftest(sa_icw_reg, cluster.vcov(sa_icw_reg, as.factor(analysis$Pole)))

#knowledge
k_fit1 <- ivreg(knowledge ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
k_se1 <- coeftest(k_fit1, cluster.vcov(k_fit1, as.factor(analysis$Pole)))

##Outcome Family Results Table## Includes all icw indices and knowledge
##Create the lists for stargazer
fit.list.icw <- list(ex_icw_reg, k_icw_reg, l_icw_reg, a_icw_reg, ha_icw_reg, sa_icw_reg, k_fit1)
se.list.icw <- list(ex_icw_reg_se[,2], k_icw_reg_se[,2], l_icw_reg_se[,2], a_icw_reg_se[,2], ha_icw_reg_se[,2], sa_icw_reg_se[,2], k_se1[,2])
p.list.icw <- list(ex_icw_reg_se[,4], k_icw_reg_se[,4], l_icw_reg_se[,4], a_icw_reg_se[,4], ha_icw_reg_se[,4], sa_icw_reg_se[,4], k_se1[,4])
depvar.list.icw <- c("Total Expenditure", "Kerosene", "Artificial Lighting", "Appliances", "Household Activity", "Satisfaction and Attitudes", "Knowledge")
omit.stats <- c("rsq","ser","f")
cov.list <- c("Legal Electrification")

##Create the table

stargazer(title = "Legal Electrification and ICW Indices",
          se = se.list.icw,
          header = F,
          fit.list.icw,
          p = p.list.icw,
          covariate.labels = cov.list,
          dep.var.labels = depvar.list.icw,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          model.names = F,
          omit = c("gender","birthplace","age","religion","caste","as.factor"),
          add.lines = list(c("Control Variables", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("N Poles", "152", "152", "152", "152", "152", "152", "152"),
                           c("N Villages", "4", "4", "4", "4", "4", "4", "4")),
          notes = "Standard errors clustered at pole level. Pole and village fixed effects included."
)


##TABLE A4##
#Total Expenditure
e_fit1 <- ivreg(total_expenditure ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
e_se1 <- coeftest(e_fit1, cluster.vcov(e_fit1, as.factor(analysis$Pole)))

#Food
e_fit2 <- ivreg(food_expenditure ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
e_se2 <- coeftest(e_fit2, cluster.vcov(e_fit2, as.factor(analysis$Pole)))

#Education
e_fit3 <- ivreg(education_expenditure ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
e_se3 <- coeftest(e_fit3, cluster.vcov(e_fit3, as.factor(analysis$Pole)))

#kerosene
e_fit4 <- ivreg(kerosene_expenditure ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
e_se4 <- coeftest(e_fit4, cluster.vcov(e_fit4, as.factor(analysis$Pole)))

#electricity
e_fit5 <- ivreg(electricity_expenditure ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
e_se5 <- coeftest(e_fit5, cluster.vcov(e_fit5, as.factor(analysis$Pole)))


##Adjust p values for multiple comparison
e_pvalues <- cbind(e_se1[,4], e_se2[,4], e_se3[,4], e_se4[,4])

#p.adjust(e_pvalues[1,], method="BH", n=length(e_pvalues[1,]))

e_pvalues_adj <- t(apply(e_pvalues, MARGIN=1, function(x) p.adjust(x, method="BH", n=length(x))))


##Create the lists for stargazer
fit.list.e <- list(e_fit1, e_fit2, e_fit3, e_fit4)
se.list.e <- list(e_se1[,2], e_se2[,2], e_se3[,2], e_se4[,2])
p.list.e <- list(e_pvalues_adj[,1], e_pvalues_adj[,2], e_pvalues_adj[,3], e_pvalues_adj[,4])
depvar.list.e <- c("Total Expenditure", "Food Expenditure", "Education Expenditure", "Kerosene Expenditure")
omit.stats <- c("rsq","ser","f")
cov.list <- c("Legal Electrification")

##Create the table
stargazer(title = "Legal Electrification and Expenditure",
          se = se.list.e,
          header = F,
          fit.list.e,
          p = p.list.e,
          covariate.labels = cov.list,
          dep.var.labels = depvar.list.e,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          model.names = F,
          omit = c("gender","birthplace","age","religion","caste","as.factor"),
          add.lines = list(c("Control Variables", "Yes", "Yes", "Yes", "Yes"),
                           c("N Poles", "152", "152", "152", "152")),
           
          notes = "Standard errors clustered at pole level. Pole fixed effects included."
)


##TABLE A5##
#lamp binary
ke_fit1 <- ivreg(kerosene_lamps ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
ke_se1 <- coeftest(ke_fit1, cluster.vcov(ke_fit1, as.factor(analysis$Pole)))

#num kerosene lamps
ke_fit2 <- ivreg(num_kerosene_lamps ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
ke_se2 <- coeftest(ke_fit2, cluster.vcov(ke_fit2, as.factor(analysis$Pole)))

#kerosene lamp hours
ke_fit3 <- ivreg(kerosene_lamp_hours ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
ke_se3 <- coeftest(ke_fit3, cluster.vcov(ke_fit3, as.factor(analysis$Pole)))

#Other use of Kerosene
ke_fit4 <- ivreg(kerosene_other ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
ke_se4 <- coeftest(ke_fit4, cluster.vcov(ke_fit4, as.factor(analysis$Pole)))


##Adjust p values for multiple comparison
ke_pvalues <- cbind(ke_se1[,4], ke_se2[,4], ke_se3[,4], ke_se4[,4])

#p.adjust(e_pvalues[1,], method="BH", n=length(e_pvalues[1,]))

ke_pvalues_adj <- t(apply(ke_pvalues, MARGIN=1, function(x) p.adjust(x, method="BH", n=length(x))))


##Create the lists for stargazer
fit.list.ke <- list(ke_fit1, ke_fit2, ke_fit3, ke_fit4)
se.list.ke <- list(ke_se1[,2], ke_se2[,2], ke_se3[,2], ke_se4[,2])
p.list.ke <- list(ke_pvalues_adj[,1], ke_pvalues_adj[,2], ke_pvalues_adj[,3], ke_pvalues_adj[,4])
depvar.list.ke <- c("Kerosene Lamp Use", "Number of Kerosene Lamps", "Kerosene Lamp Hours", "Other Use of Kerosene")
omit.stats <- c("rsq","ser","f")
cov.list <- c("Legal Electrification")

##Create the table
stargazer(title = "Legal Electrification and Kerosene Use",
          se = se.list.ke,
          header = F,
          fit.list.ke,
          p = p.list.ke,
          covariate.labels = cov.list,
          dep.var.labels = depvar.list.ke,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          model.names = F,
          omit = c("gender","birthplace","age","religion","caste","as.factor"),
          add.lines = list(c("Control Variables", "Yes", "Yes", "Yes", "Yes"),
                           c("N Poles", "152", "152", "152", "152")),
          
          notes = "Standard errors clustered at pole level. Pole fixed effects included."
)


##TABLE A6##
#Lighting hours
l_fit1 <- ivreg(lighting_hours ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
l_se1 <- coeftest(l_fit1, cluster.vcov(l_fit1, as.factor(analysis$Pole)))

#child lighting hours
l_fit2 <- ivreg(child_lighting ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
l_se2 <- coeftest(l_fit2, cluster.vcov(l_fit2, as.factor(analysis$Pole)))

#adult lighting hours
l_fit3 <- ivreg(adult_lighting ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
l_se3 <- coeftest(l_fit3, cluster.vcov(l_fit3, as.factor(analysis$Pole)))


##Adjust p values for multiple comparison
l_pvalues <- cbind(l_se1[,4], l_se2[,4], l_se3[,4])

#p.adjust(e_pvalues[1,], method="BH", n=length(e_pvalues[1,]))

l_pvalues_adj <- t(apply(l_pvalues, MARGIN=1, function(x) p.adjust(x, method="BH", n=length(x))))

##Create the lists for stargazer
fit.list.l <- list(l_fit1, l_fit2, l_fit3)
se.list.l <- list(l_se1[,2], l_se2[,2], l_se3[,2])
p.list.l <- list(l_pvalues_adj[,1], l_pvalues_adj[,2], l_pvalues_adj[,3])
depvar.list.l <- c("Household Light", "Child Light", "Adult Light")
omit.stats <- c("rsq","ser","f")
cov.list <- c("Legal Electrification")

##Create the table
stargazer(title = "Legal Electrification and Lighting Use",
          se = se.list.l,
          header = F,
          fit.list.l,
          p = p.list.l,
          covariate.labels = cov.list,
          dep.var.labels = depvar.list.l,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          model.names = F,
          omit = c("gender","birthplace","age","religion","caste","as.factor"),
          add.lines = list(c("Control Variables", "Yes", "Yes", "Yes"),
                           c("N Poles", "152", "152", "152")),
           
          notes = "Standard errors clustered at pole level. Pole fixed effects included."
)


##TABLE A7##
#adult activity
ha_fit1 <- ivreg(adult_activity ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
ha_se1 <- coeftest(ha_fit1, cluster.vcov(ha_fit1, as.factor(analysis$Pole)))

#child activity
ha_fit2 <- ivreg(child_activity ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
ha_se2 <- coeftest(ha_fit2, cluster.vcov(ha_fit2, as.factor(analysis$Pole)))

##Adjust p values for multiple comparison
ha_pvalues <- cbind(ha_se1[,4], ha_se2[,4])

#p.adjust(e_pvalues[1,], method="BH", n=length(e_pvalues[1,]))

ha_pvalues_adj <- t(apply(ha_pvalues, MARGIN=1, function(x) p.adjust(x, method="BH", n=length(x))))

##Create the lists for stargazer
fit.list.ha <- list(ha_fit1, ha_fit2)
se.list.ha <- list(ha_se1[,2], ha_se2[,2])
p.list.ha <- list(ha_pvalues_adj[,1], ha_pvalues_adj[,2])
depvar.list.ha <- c("Adult Activity", "Child Activity")
omit.stats <- c("rsq","ser","f")
cov.list <- c("Legal Electrification")

##Create the table
stargazer(title = "Legal Electrification and Household Activity",
          se = se.list.ha,
          header = F,
          fit.list.ha,
          p = p.list.ha,
          covariate.labels = cov.list,
          dep.var.labels = depvar.list.ha,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          model.names = F,
          omit = c("gender","birthplace","age","religion","caste","as.factor"),
          add.lines = list(c("Control Variables", "Yes", "Yes"),
                           c("N Poles", "152", "152")),
          notes = "Standard errors clustered at pole level. Pole fixed effects included."
)

##TABLE A8##
a_fit1 <- ivreg(appliances ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
a_se1 <- coeftest(a_fit1, cluster.vcov(a_fit1, as.factor(analysis$Pole)))

#appliance use
a_fit2 <- ivreg(appliance_use ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
a_se2 <- coeftest(a_fit2, cluster.vcov(a_fit2, as.factor(analysis$Pole)))

##Adjust p values for multiple comparison
a_pvalues <- cbind(a_se1[,4], a_se2[,4])

#p.adjust(e_pvalues[1,], method="BH", n=length(e_pvalues[1,]))

a_pvalues_adj <- t(apply(a_pvalues, MARGIN=1, function(x) p.adjust(x, method="BH", n=length(x))))

##Create the lists for stargazer
fit.list.a <- list(a_fit1, a_fit2)
se.list.a <- list(a_se1[,2], a_se2[,2])
p.list.a <- list(a_pvalues_adj[,1], a_pvalues_adj[,2])
depvar.list.a <- c("Number of Appliances", "Appliance Use")
omit.stats <- c("rsq","ser","f")
cov.list <- c("Legal Electrification")

##Create the table
stargazer(title = "Legal Electrification and Appliances",
          se = se.list.a,
          header = F,
          fit.list.a,
          p = p.list.a,
          covariate.labels = cov.list,
          dep.var.labels = depvar.list.a,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          model.names = F,
          omit = c("gender","birthplace","age","religion","caste","as.factor"),
          add.lines = list(c("Control Variables", "Yes", "Yes"),
                           c("N Poles", "152", "152")),
           
          notes = "Standard errors clustered at pole level. Pole fixed effects included."
)


##TABLE A9##
#reliability
sa_fit1 <- ivreg(satisfaction_reliability ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
sa_se1 <- coeftest(sa_fit1, cluster.vcov(sa_fit1, as.factor(analysis$Pole)))

#cost
sa_fit2 <- ivreg(satisfaction_cost ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
sa_se2 <- coeftest(sa_fit2, cluster.vcov(sa_fit2, as.factor(analysis$Pole)))

#safety
sa_fit3 <- ivreg(satisfaction_safety ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
sa_se3 <- coeftest(sa_fit3, cluster.vcov(sa_fit3, as.factor(analysis$Pole)))

#brightness
sa_fit4 <- ivreg(satisfaction_brightness ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
sa_se4 <- coeftest(sa_fit4, cluster.vcov(sa_fit4, as.factor(analysis$Pole)))


#overall satisfaction
sa_fit5 <- ivreg(satisfaction ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
sa_se5 <- coeftest(sa_fit5, cluster.vcov(sa_fit5, as.factor(analysis$Pole)))

#change in satisfaction
sa_fit6 <- ivreg(satisfaction_chng ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
sa_se6 <- coeftest(sa_fit6, cluster.vcov(sa_fit6, as.factor(analysis$Pole)))

#electricity value
sa_fit7 <- ivreg(elec_value ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
sa_se7 <- coeftest(sa_fit7, cluster.vcov(sa_fit7, as.factor(analysis$Pole)))

#satisfaction with business
sa_fit8 <- ivreg(satisfaction_business ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
sa_se8 <- coeftest(sa_fit8, cluster.vcov(sa_fit8, as.factor(analysis$Pole)))

#increase in business income
sa_fit9 <- ivreg(income_increase ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
sa_se9 <- coeftest(sa_fit9, cluster.vcov(sa_fit9, as.factor(analysis$Pole)))

#interest in business
sa_fit10 <- ivreg(business_interest ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
sa_se10 <- coeftest(sa_fit10, cluster.vcov(sa_fit10, as.factor(analysis$Pole)))

#aspirations/modernity
sa_fit11 <- ivreg(aspirations ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
sa_se11 <- coeftest(sa_fit11, cluster.vcov(sa_fit11, as.factor(analysis$Pole)))


##Adjust P values
sa_pvalues <- cbind(sa_se1[,4], sa_se2[,4], sa_se3[,4], sa_se4[,4], sa_se5[,4], sa_se6[,4], sa_se7[,4], sa_se8[,4], sa_se9[,4], 
                    sa_se10[,4], sa_se11[,4])

#p.adjust(e_pvalues[1,], method="BH", n=length(e_pvalues[1,]))

sa_pvalues_adj <- t(apply(sa_pvalues, MARGIN=1, function(x) p.adjust(x, method="BH", n=length(x))))




##Create the lists for stargazer


fit.list.sa1 <- list(sa_fit1, sa_fit2, sa_fit3, sa_fit4)
se.list.sa1 <- list(sa_se1[,2], sa_se2[,2], sa_se3[,2], sa_se4[,2])
p.list.sa1 <- list(sa_pvalues_adj[,1], sa_pvalues_adj[,2], sa_pvalues_adj[,3], sa_pvalues_adj[,4])
depvar.list.sa1 <- c("Reliability", "Cost", "Safety", "Brightness")
cov.list <- c("Legal Electrification")

#Create the table
stargazer(title = "Legal Electrification and Satisfaction and Attitudes",
          se = se.list.sa1,
          header = F,
          fit.list.sa1,
          p = p.list.sa1,
          covariate.labels = cov.list,
          dep.var.labels = depvar.list.sa1,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          model.names = F,
          omit = c("gender","birthplace","age","religion","caste","as.factor"),
          add.lines = list(c("Control Variables", "Yes", "Yes", "Yes", "Yes"),
                           c("N Poles", "152", "152", "152", "152")),
          notes = "Standard errors clustered at pole level. Pole fixed effects included."
)


##TABLE A10##
fit.list.sa2 <- list(sa_fit5, sa_fit6, sa_fit7)
se.list.sa2 <- list(sa_se5[,2], sa_se6[,2], sa_se7[,2])
p.list.sa2 <- list(sa_pvalues_adj[,5], sa_pvalues_adj[,6], sa_pvalues_adj[,7])
depvar.list.sa2 <- c("Overall Satisfaction", "Change in Satisfaction", "Electricity Value")

#Create the table
stargazer(title = "Legal Electrification and Satisfaction and Attitudes",
          se = se.list.sa2,
          header = F,
          fit.list.sa2,
          p = p.list.sa2,
          covariate.labels = cov.list,
          dep.var.labels = depvar.list.sa2,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          model.names = F,
          omit = c("gender","birthplace","age","religion","caste","as.factor"),
          add.lines = list(c("Control Variables", "Yes", "Yes", "Yes"),
                           c("N Poles", "152", "152", "152")),
          notes = "Standard errors clustered at pole level. Pole fixed effects included."
)


##TABLE A11##
fit.list.sa3 <- list(sa_fit9, sa_fit10, sa_fit11)
se.list.sa3 <- list(sa_se9[,2], sa_se10[,2], sa_se11[,2])
p.list.sa3 <- list(sa_pvalues_adj[,9], sa_pvalues_adj[,10], sa_pvalues_adj[,11])
depvar.list.sa3 <- c("Increase in Income", "Business Aspiration", "General Aspirations")

#Create the table
stargazer(title = "Legal Electrification and Satisfaction and Attitudes",
          se = se.list.sa3,
          header = F,
          fit.list.sa3,
          p = p.list.sa3,
          covariate.labels = cov.list,
          dep.var.labels = depvar.list.sa3,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          model.names = F,
          omit = c("gender","birthplace","age","religion","caste","as.factor"),
          add.lines = list(c("Control Variables", "Yes", "Yes", "Yes"),
                           c("N Poles", "152", "152", "152")),
           
          notes = "Standard errors clustered at pole level. Pole fixed effects included."
)

##TABLE A12##
k_fit1 <- ivreg(knowledge ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
k_se1 <- coeftest(k_fit1, cluster.vcov(k_fit1, as.factor(analysis$Pole)))

#no need for padjust- only 1 model

##Create the lists for stargazer
fit.list.k <- list(k_fit1)
se.list.k <- list(k_se1[,2])
omit.stats <- c("rsq","ser","f")
cov.list <- c("Legal Electrification")


##Create the table
stargazer(title = "Legal Electrification and Knowledge",
          se = se.list.k,
          header = F,
          fit.list.k,
          covariate.labels = cov.list,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          model.names = F,
          omit = c("gender","birthplace","age","religion","caste","as.factor"),
          add.lines = list(c("Control Variables", "Yes"),
                           c("N Poles", "152")),
          dep.var.labels = "Knowledge",
          notes = "Standard errors clustered at pole level. Pole fixed effects included."
)


##TABLE A13##
#Total Expenditure
e_fit1_rf <- lm(total_expenditure ~ forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
e_se1_rf <- coeftest(e_fit1_rf, cluster.vcov(e_fit1_rf, as.factor(analysis$Pole)))

#Food
e_fit2_rf <- lm(food_expenditure ~ forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
e_se2_rf <- coeftest(e_fit2_rf, cluster.vcov(e_fit2_rf, as.factor(analysis$Pole)))

#Education
e_fit3_rf <- lm(education_expenditure ~ forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
e_se3_rf <- coeftest(e_fit3_rf, cluster.vcov(e_fit3_rf, as.factor(analysis$Pole)))

#kerosene
e_fit4_rf <- lm(kerosene_expenditure ~ forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
e_se4_rf <- coeftest(e_fit4_rf, cluster.vcov(e_fit4_rf, as.factor(analysis$Pole)))




##Adjust p values for multiple comparison
e_pvalues_rf <- cbind(e_se1_rf[,4], e_se2_rf[,4], e_se3_rf[,4], e_se4_rf[,4])

#p.adjust(e_pvalues[1,], method="BH", n=length(e_pvalues[1,]))

e_pvalues_adj_rf <- t(apply(e_pvalues_rf, MARGIN=1, function(x) p.adjust(x, method="BH", n=length(x))))


##Create the lists for stargazer
fit.list.e_rf <- list(e_fit1_rf, e_fit2_rf, e_fit3_rf, e_fit4_rf)
se.list.e_rf <- list(e_se1_rf[,2], e_se2_rf[,2], e_se3_rf[,2], e_se4_rf[,2])
p.list.e_rf <- list(e_pvalues_adj_rf[,1], e_pvalues_adj_rf[,2], e_pvalues_adj_rf[,3], e_pvalues_adj_rf[,4])
depvar.list.e <- c("Total Expenditure", "Food Expenditure", "Education Expenditure", "Kerosene Expenditure")
omit.stats <- c("rsq","ser","f")
cov.list.rf <- c("Distance (Binary)")

##Create the table
stargazer(title = "Reduced Form Analysis of Legal Electrification and Expenditure",
          se = se.list.e_rf,
          header = F,
          fit.list.e_rf,
          p = p.list.e_rf,
          covariate.labels = cov.list.rf,
          dep.var.labels = depvar.list.e,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          model.names = F,
          omit = c("gender","birthplace","age","religion","caste","as.factor"),
          add.lines = list(c("Control Variables", "Yes", "Yes", "Yes", "Yes"),
                           c("N Poles", "152", "152", "152", "152")),
           
          notes = "Standard errors clustered at pole level. Pole fixed effects included."
)

##TABLE A14##
#lamp binary
ke_fit1_rf <- lm(kerosene_lamps ~ forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
ke_se1_rf <- coeftest(ke_fit1_rf, cluster.vcov(ke_fit1_rf, as.factor(analysis$Pole)))

#num kerosene lamps
ke_fit2_rf <- lm(num_kerosene_lamps ~  forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
ke_se2_rf <- coeftest(ke_fit2_rf, cluster.vcov(ke_fit2_rf, as.factor(analysis$Pole)))

#kerosene lamp hours
ke_fit3_rf <- lm(kerosene_lamp_hours ~  forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
ke_se3_rf <- coeftest(ke_fit3_rf, cluster.vcov(ke_fit3_rf, as.factor(analysis$Pole)))

#Other use of Kerosene
ke_fit4_rf <- lm(kerosene_other ~  forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
ke_se4_rf <- coeftest(ke_fit4_rf, cluster.vcov(ke_fit4_rf, as.factor(analysis$Pole)))


##Adjust p values for multiple comparison
ke_pvalues_rf <- cbind(ke_se1_rf[,4], ke_se2_rf[,4], ke_se3_rf[,4], ke_se4_rf[,4])

#p.adjust(e_pvalues[1,], method="BH", n=length(e_pvalues[1,]))

ke_pvalues_adj_rf <- t(apply(ke_pvalues_rf, MARGIN=1, function(x) p.adjust(x, method="BH", n=length(x))))


##Create the lists for stargazer
fit.list.ke_rf <- list(ke_fit1_rf, ke_fit2_rf, ke_fit3_rf, ke_fit4_rf)
se.list.ke_rf <- list(ke_se1_rf[,2], ke_se2_rf[,2], ke_se3_rf[,2], ke_se4_rf[,2])
p.list.ke_rf <- list(ke_pvalues_adj_rf[,1], ke_pvalues_adj_rf[,2], ke_pvalues_adj_rf[,3], ke_pvalues_adj_rf[,4])
depvar.list.ke <- c("Kerosene Lamp Use", "Number of Kerosene Lamps", "Kerosene Lamp Hours", "Other Use of Kerosene")
omit.stats <- c("rsq","ser","f")

##Create the table
stargazer(title = "Reduced Form Legal Electrification and Kerosene Use",
          se = se.list.ke_rf,
          header = F,
          fit.list.ke_rf,
          p = p.list.ke_rf,
          covariate.labels = cov.list.rf,
          dep.var.labels = depvar.list.ke,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          model.names = F,
          omit = c("gender","birthplace","age","religion","caste","as.factor"),
          add.lines = list(c("Control Variables", "Yes", "Yes", "Yes", "Yes"),
                           c("N Poles", "152", "152", "152", "152")),
          
          notes = "Standard errors clustered at pole level. Pole fixed effects included."
)

##TABLE A15##
#Lighting hours
l_fit1_rf <- lm(lighting_hours ~  forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
l_se1_rf <- coeftest(l_fit1_rf, cluster.vcov(l_fit1_rf, as.factor(analysis$Pole)))

#child lighting hours
l_fit2_rf <- lm(child_lighting ~  forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
l_se2_rf <- coeftest(l_fit2_rf, cluster.vcov(l_fit2_rf, as.factor(analysis$Pole)))

#adult lighting hours
l_fit3_rf <- lm(adult_lighting ~  forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
l_se3_rf <- coeftest(l_fit3_rf, cluster.vcov(l_fit3_rf, as.factor(analysis$Pole)))


##Adjust p values for multiple comparison
l_pvalues_rf <- cbind(l_se1_rf[,4], l_se2_rf[,4], l_se3_rf[,4])

#p.adjust(e_pvalues[1,], method="BH", n=length(e_pvalues[1,]))

l_pvalues_adj_rf <- t(apply(l_pvalues_rf, MARGIN=1, function(x) p.adjust(x, method="BH", n=length(x))))

##Create the lists for stargazer
fit.list.l_rf <- list(l_fit1_rf, l_fit2_rf, l_fit3_rf)
se.list.l_rf <- list(l_se1_rf[,2], l_se2_rf[,2], l_se3_rf[,2])
p.list.l_rf <- list(l_pvalues_adj_rf[,1], l_pvalues_adj_rf[,2], l_pvalues_adj_rf[,3])
depvar.list.l <- c("Household Light", "Child Light", "Adult Light")
omit.stats <- c("rsq","ser","f")

##Create the table
stargazer(title = "Reduced Form Legal Electrification and Lighting Use",
          se = se.list.l_rf,
          header = F,
          fit.list.l_rf,
          p = p.list.l_rf,
          covariate.labels = cov.list.rf,
          dep.var.labels = depvar.list.l,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          model.names = F,
          omit = c("gender","birthplace","age","religion","caste","as.factor"),
          add.lines = list(c("Control Variables", "Yes", "Yes", "Yes"),
                           c("N Poles", "152", "152", "152")),
           
          notes = "Standard errors clustered at pole level. Pole fixed effects included."
)

##TABLE A16##
#adult activity
ha_fit1_rf <- lm(adult_activity ~  forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
ha_se1_rf <- coeftest(ha_fit1_rf, cluster.vcov(ha_fit1_rf, as.factor(analysis$Pole)))

#child activity
ha_fit2_rf <- lm(child_activity ~  forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
ha_se2_rf <- coeftest(ha_fit2_rf, cluster.vcov(ha_fit2_rf, as.factor(analysis$Pole)))

##Adjust p values for multiple comparison
ha_pvalues_rf <- cbind(ha_se1_rf[,4], ha_se2_rf[,4])

#p.adjust(e_pvalues[1,], method="BH", n=length(e_pvalues[1,]))

ha_pvalues_adj_rf <- t(apply(ha_pvalues_rf, MARGIN=1, function(x) p.adjust(x, method="BH", n=length(x))))

##Create the lists for stargazer
fit.list.ha_rf <- list(ha_fit1_rf, ha_fit2_rf)
se.list.ha_rf <- list(ha_se1_rf[,2], ha_se2_rf[,2])
p.list.ha_rf <- list(ha_pvalues_adj_rf[,1], ha_pvalues_adj_rf[,2])
depvar.list.ha <- c("Adult Activity", "Child Activity")
omit.stats <- c("rsq","ser","f")

##Create the table
stargazer(title = "Reduced Form Legal Electrification and Household Activity",
          se = se.list.ha_rf,
          header = F,
          fit.list.ha_rf,
          p = p.list.ha_rf,
          covariate.labels = cov.list.rf,
          dep.var.labels = depvar.list.ha,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          model.names = F,
          omit = c("gender","birthplace","age","religion","caste","as.factor"),
          add.lines = list(c("Control Variables", "Yes", "Yes"),
                           c("N Poles", "152", "152")),
          notes = "Standard errors clustered at pole level. Pole fixed effects included."
)

##TABLE A17##
##Appliances##
a_fit1_rf <- lm(appliances ~  forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
a_se1_rf <- coeftest(a_fit1_rf, cluster.vcov(a_fit1_rf, as.factor(analysis$Pole)))

#appliance use
a_fit2_rf <- lm(appliance_use ~  forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
a_se2_rf <- coeftest(a_fit2_rf, cluster.vcov(a_fit2_rf, as.factor(analysis$Pole)))

##Adjust p values for multiple comparison
a_pvalues_rf <- cbind(a_se1_rf[,4], a_se2_rf[,4])

#p.adjust(e_pvalues[1,], method="BH", n=length(e_pvalues[1,]))

a_pvalues_adj_rf <- t(apply(a_pvalues_rf, MARGIN=1, function(x) p.adjust(x, method="BH", n=length(x))))

##Create the lists for stargazer
fit.list.a_rf <- list(a_fit1_rf, a_fit2_rf)
se.list.a_rf <- list(a_se1_rf[,2], a_se2_rf[,2])
p.list.a_rf <- list(a_pvalues_adj_rf[,1], a_pvalues_adj_rf[,2])
depvar.list.a <- c("Number of Appliances", "Appliance Use")
omit.stats <- c("rsq","ser","f")

##Create the table
stargazer(title = "Reduced Form Legal Electrification and Appliances",
          se = se.list.a_rf,
          header = F,
          fit.list.a_rf,
          p = p.list.a_rf,
          covariate.labels = cov.list.rf,
          dep.var.labels = depvar.list.a,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          model.names = F,
          omit = c("gender","birthplace","age","religion","caste","as.factor"),
          add.lines = list(c("Control Variables", "Yes", "Yes"),
                           c("N Poles", "152", "152")),
           
          notes = "Standard errors clustered at pole level. Pole fixed effects included."
)


##TABLE A18##
##Satisfaction and Attitudes##
#reliability
sa_fit1_rf <- lm(satisfaction_reliability ~  forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
sa_se1_rf <- coeftest(sa_fit1_rf, cluster.vcov(sa_fit1_rf, as.factor(analysis$Pole)))

#cost
sa_fit2_rf <- lm(satisfaction_cost ~  forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
sa_se2_rf <- coeftest(sa_fit2_rf, cluster.vcov(sa_fit2_rf, as.factor(analysis$Pole)))

#safety
sa_fit3_rf <- lm(satisfaction_safety ~  forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
sa_se3_rf <- coeftest(sa_fit3_rf, cluster.vcov(sa_fit3_rf, as.factor(analysis$Pole)))

#brightness
sa_fit4_rf <- lm(satisfaction_brightness ~  forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
sa_se4_rf <- coeftest(sa_fit4_rf, cluster.vcov(sa_fit4_rf, as.factor(analysis$Pole)))


#overall satisfaction
sa_fit5_rf <- lm(satisfaction ~  forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
sa_se5_rf <- coeftest(sa_fit5_rf, cluster.vcov(sa_fit5_rf, as.factor(analysis$Pole)))

#change in satisfaction
sa_fit6_rf <- lm(satisfaction_chng ~  forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
sa_se6_rf <- coeftest(sa_fit6_rf, cluster.vcov(sa_fit6_rf, as.factor(analysis$Pole)))

#electricity value
sa_fit7_rf <- lm(elec_value ~  forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
sa_se7_rf <- coeftest(sa_fit7_rf, cluster.vcov(sa_fit7_rf, as.factor(analysis$Pole)))

#satisfaction with business
sa_fit8_rf <- lm(satisfaction_business ~  forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
sa_se8_rf <- coeftest(sa_fit8_rf, cluster.vcov(sa_fit8_rf, as.factor(analysis$Pole)))

#increase in business income
sa_fit9_rf <- lm(income_increase ~  forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
sa_se9_rf <- coeftest(sa_fit9_rf, cluster.vcov(sa_fit9_rf, as.factor(analysis$Pole)))

#interest in business
sa_fit10_rf <- lm(business_interest ~  forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
sa_se10_rf <- coeftest(sa_fit10_rf, cluster.vcov(sa_fit10_rf, as.factor(analysis$Pole)))

#aspirations/modernity
sa_fit11_rf <- lm(aspirations ~  forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
sa_se11_rf <- coeftest(sa_fit11_rf, cluster.vcov(sa_fit11_rf, as.factor(analysis$Pole)))


##Adjust P values
sa_pvalues_rf <- cbind(sa_se1_rf[,4], sa_se2_rf[,4], sa_se3_rf[,4], sa_se4_rf[,4], sa_se5_rf[,4], sa_se6_rf[,4], sa_se7_rf[,4], sa_se8_rf[,4], sa_se9_rf[,4], 
                       sa_se10_rf[,4], sa_se11_rf[,4])

#p.adjust(e_pvalues[1,], method="BH", n=length(e_pvalues[1,]))

sa_pvalues_adj_rf <- t(apply(sa_pvalues_rf, MARGIN=1, function(x) p.adjust(x, method="BH", n=length(x))))




##Create the lists for stargazer

#Table sa1
fit.list.sa1_rf <- list(sa_fit1_rf, sa_fit2_rf, sa_fit3_rf, sa_fit4_rf)
se.list.sa1_rf <- list(sa_se1_rf[,2], sa_se2_rf[,2], sa_se3_rf[,2], sa_se4_rf[,2])
p.list.sa1_rf <- list(sa_pvalues_adj_rf[,1], sa_pvalues_adj_rf[,2], sa_pvalues_adj_rf[,3], sa_pvalues_adj_rf[,4])
depvar.list.sa1 <- c("Reliability", "Cost", "Safety", "Brightness")

#Create the table
stargazer(title = "Reduced Form Legal Electrification and Satisfaction and Attitudes",
          se = se.list.sa1_rf,
          header = F,
          fit.list.sa1_rf,
          p = p.list.sa1_rf,
          covariate.labels = cov.list.rf,
          dep.var.labels = depvar.list.sa1,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          model.names = F,
          omit = c("gender","birthplace","age","religion","caste","as.factor"),
          add.lines = list(c("Control Variables", "Yes", "Yes", "Yes", "Yes"),
                           c("N Poles", "152", "152", "152", "152")),
           
          notes = "Standard errors clustered at pole level. Pole fixed effects included."
)


##TABLE A19##
fit.list.sa2_rf <- list(sa_fit5_rf, sa_fit6_rf, sa_fit7_rf)
se.list.sa2_rf <- list(sa_se5_rf[,2], sa_se6_rf[,2], sa_se7_rf[,2])
p.list.sa2_rf <- list(sa_pvalues_adj_rf[,5], sa_pvalues_adj_rf[,6], sa_pvalues_adj_rf[,7])
depvar.list.sa2 <- c("Overall Satisfaction", "Change in Satisfaction", "Electricity Value")

#Create the table
stargazer(title = "Reduced Form Legal Electrification and Satisfaction and Attitudes",
          se = se.list.sa2_rf,
          header = F,
          fit.list.sa2_rf,
          p = p.list.sa2_rf,
          covariate.labels = cov.list.rf,
          dep.var.labels = depvar.list.sa2,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          model.names = F,
          omit = c("gender","birthplace","age","religion","caste","as.factor"),
          add.lines = list(c("Control Variables", "Yes", "Yes", "Yes"),
                           c("N Poles", "152", "152", "152")),
           
          notes = "Standard errors clustered at pole level. Pole fixed effects included."
)



##TABLE A20##
fit.list.sa3_rf <- list(sa_fit9_rf, sa_fit10_rf, sa_fit11_rf)
se.list.sa3_rf <- list(sa_se9_rf[,2], sa_se10_rf[,2], sa_se11_rf[,2])
p.list.sa3_rf <- list(sa_pvalues_adj_rf[,9], sa_pvalues_adj_rf[,10], sa_pvalues_adj_rf[,11])
depvar.list.sa3 <- c("Increase in Income", "Business Aspiration", "General Aspirations")

#Create the table
stargazer(title = "Reduced Form Legal Electrification and Satisfaction and Attitudes",
          se = se.list.sa3_rf,
          header = F,
          fit.list.sa3_rf,
          p = p.list.sa3_rf,
          covariate.labels = cov.list.rf,
          dep.var.labels = depvar.list.sa3,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          model.names = F,
          omit = c("gender","birthplace","age","religion","caste","as.factor"),
          add.lines = list(c("Control Variables", "Yes", "Yes", "Yes"),
                           c("N Poles", "152", "152", "152")),
           
          notes = "Standard errors clustered at pole level. Pole fixed effects included."
)



##TABLE A21##
k_fit1_rf <- lm(knowledge ~  forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
k_se1_rf <- coeftest(k_fit1_rf, cluster.vcov(k_fit1_rf, as.factor(analysis$Pole)))

##Create the lists for stargazer
fit.list.k_rf <- list(k_fit1_rf)
se.list.k_rf <- list(k_se1_rf[,2])
omit.stats <- c("rsq","ser","f")



##Create the table
stargazer(title = "Legal Electrification and Knowledge",
          se = se.list.k_rf,
          header = F,
          fit.list.k_rf,
          covariate.labels = cov.list.rf,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          model.names = F,
          omit = c("gender","birthplace","age","religion","caste","as.factor"),
          add.lines = list(c("Control Variables", "Yes"),
                           c("N Poles", "152")),
          dep.var.labels = "Knowledge",
          notes = "Standard errors clustered at pole level. Pole fixed effects included."
)


##TABLE A22##
e_fit1v <- ivreg(total_expenditure ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) +as.factor(village) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole) +as.factor(village), data=analysis)
e_se1v <- coeftest(e_fit1v, cluster.vcov(e_fit1v, as.factor(analysis$Pole)))

#Food
e_fit2v <- ivreg(food_expenditure ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) +as.factor(village) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole) +as.factor(village), data=analysis)
e_se2v <- coeftest(e_fit2v, cluster.vcov(e_fit2v, as.factor(analysis$Pole)))

#Education
e_fit3v <- ivreg(education_expenditure ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village), data=analysis)
e_se3v <- coeftest(e_fit3v, cluster.vcov(e_fit3v, as.factor(analysis$Pole)))

#kerosene
e_fit4v <- ivreg(kerosene_expenditure ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village), data=analysis)
e_se4v <- coeftest(e_fit4v, cluster.vcov(e_fit4v, as.factor(analysis$Pole)))




##Adjust p values for multiple comparison
e_pvaluesv <- cbind(e_se1v[,4], e_se2v[,4], e_se3v[,4], e_se4v[,4])


e_pvalues_adjv <- t(apply(e_pvaluesv, MARGIN=1, function(x) p.adjust(x, method="BH", n=length(x))))


##Create the lists for stargazer
fit.list.ev <- list(e_fit1v, e_fit2v, e_fit3v, e_fit4v)
se.list.ev <- list(e_se1v[,2], e_se2v[,2], e_se3v[,2], e_se4v[,2])
p.list.ev <- list(e_pvalues_adjv[,1], e_pvalues_adjv[,2], e_pvalues_adjv[,3], e_pvalues_adjv[,4])
depvar.list.e <- c("Total Expenditure", "Food Expenditure", "Education Expenditure", "Kerosene Expenditure")
omit.stats <- c("rsq","ser","f")
cov.list <- c("Legal Electrification")

##Create the table
stargazer(title = "Legal Electrification and Expenditure",
          se = se.list.ev,
          header = F,
          fit.list.ev,
          p = p.list.ev,
          covariate.labels = cov.list,
          dep.var.labels = depvar.list.e,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          model.names = F,
          omit = c("gender","birthplace","age","religion","caste","as.factor"),
          add.lines = list(c("Control Variables", "Yes", "Yes", "Yes", "Yes"),
                           c("N Poles", "152", "152", "152", "152"),
                           c("N Villages", "4", "4", "4", "4")),
           
          notes = "Standard errors clustered at pole level. Pole and village fixed effects included."
)


##TABLE A23##
ke_fit1v <- ivreg(kerosene_lamps ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village), data=analysis)
ke_se1v <- coeftest(ke_fit1v, cluster.vcov(ke_fit1v, as.factor(analysis$Pole)))

#num kerosene lamps
ke_fit2v <- ivreg(num_kerosene_lamps ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village)| forcing + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village), data=analysis)
ke_se2v <- coeftest(ke_fit2v, cluster.vcov(ke_fit2v, as.factor(analysis$Pole)))

#kerosene lamp hours
ke_fit3v <- ivreg(kerosene_lamp_hours ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village), data=analysis)
ke_se3v <- coeftest(ke_fit3v, cluster.vcov(ke_fit3v, as.factor(analysis$Pole)))

#Other use of Kerosene
ke_fit4v <- ivreg(kerosene_other ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village), data=analysis)
ke_se4v <- coeftest(ke_fit4v, cluster.vcov(ke_fit4v, as.factor(analysis$Pole)))


##Adjust p values for multiple comparison
ke_pvaluesv <- cbind(ke_se1v[,4], ke_se2v[,4], ke_se3v[,4], ke_se4v[,4])

#p.adjust(e_pvalues[1,], method="BH", n=length(e_pvalues[1,]))

ke_pvalues_adjv <- t(apply(ke_pvaluesv, MARGIN=1, function(x) p.adjust(x, method="BH", n=length(x))))


##Create the lists for stargazer
fit.list.kev <- list(ke_fit1v, ke_fit2v, ke_fit3v, ke_fit4v)
se.list.kev <- list(ke_se1v[,2], ke_se2v[,2], ke_se3v[,2], ke_se4v[,2])
p.list.kev <- list(ke_pvalues_adjv[,1], ke_pvalues_adjv[,2], ke_pvalues_adjv[,3], ke_pvalues_adjv[,4])
depvar.list.ke <- c("Kerosene Lamp Use", "Number of Kerosene Lamps", "Kerosene Lamp Hours", "Other Use of Kerosene")
omit.stats <- c("rsq","ser","f")
cov.list <- c("Legal Electrification")

##Create the table
stargazer(title = "Legal Electrification and Kerosene Use",
          se = se.list.kev,
          header = F,
          fit.list.kev,
          p = p.list.kev,
          covariate.labels = cov.list,
          dep.var.labels = depvar.list.ke,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          model.names = F,
          omit = c("gender","birthplace","age","religion","caste","as.factor"),
          add.lines = list(c("Control Variables", "Yes", "Yes", "Yes", "Yes"),
                           c("N Poles", "152", "152", "152", "152"),
                           c("N Villages", "4", "4", "4", "4")),
          
          notes = "Standard errors clustered at pole level. Pole and village fixed effects included."
)


##TABLE A24##
#Lighting hours
l_fit1v <- ivreg(lighting_hours ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village), data=analysis)
l_se1v <- coeftest(l_fit1v, cluster.vcov(l_fit1v, as.factor(analysis$Pole)))

#child lighting hours
l_fit2v <- ivreg(child_lighting ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village), data=analysis)
l_se2v <- coeftest(l_fit2v, cluster.vcov(l_fit2v, as.factor(analysis$Pole)))

#adult lighting hours
l_fit3v <- ivreg(adult_lighting ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village), data=analysis)
l_se3v <- coeftest(l_fit3v, cluster.vcov(l_fit3v, as.factor(analysis$Pole)))


##Adjust p values for multiple comparison
l_pvaluesv <- cbind(l_se1v[,4], l_se2v[,4], l_se3v[,4])


l_pvalues_adjv <- t(apply(l_pvaluesv, MARGIN=1, function(x) p.adjust(x, method="BH", n=length(x))))

##Create the lists for stargazer
fit.list.lv <- list(l_fit1v, l_fit2v, l_fit3v)
se.list.lv <- list(l_se1v[,2], l_se2v[,2], l_se3v[,2])
p.list.lv <- list(l_pvalues_adjv[,1], l_pvalues_adjv[,2], l_pvalues_adjv[,3])
depvar.list.l <- c("Household Light", "Child Light", "Adult Light")
omit.stats <- c("rsq","ser","f")
cov.list <- c("Legal Electrification")

##Create the table
stargazer(title = "Legal Electrification and Lighting Use",
          se = se.list.lv,
          header = F,
          fit.list.lv,
          p = p.list.lv,
          covariate.labels = cov.list,
          dep.var.labels = depvar.list.l,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          model.names = F,
          omit = c("gender","birthplace","age","religion","caste","as.factor"),
          add.lines = list(c("Control Variables", "Yes", "Yes", "Yes"),
                           c("N Poles", "152", "152", "152"),
                           c("N Villages", "4", "4", "4")),
           
          notes = "Standard errors clustered at pole level. Pole and village fixed effects included."
)

##TABLE A25##
#adult activity
ha_fit1v <- ivreg(adult_activity ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village), data=analysis)
ha_se1v <- coeftest(ha_fit1v, cluster.vcov(ha_fit1v, as.factor(analysis$Pole)))

#child activity
ha_fit2v <- ivreg(child_activity ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village), data=analysis)
ha_se2v <- coeftest(ha_fit2v, cluster.vcov(ha_fit2v, as.factor(analysis$Pole)))

##Adjust p values for multiple comparison
ha_pvaluesv <- cbind(ha_se1v[,4], ha_se2v[,4])


ha_pvalues_adjv <- t(apply(ha_pvaluesv, MARGIN=1, function(x) p.adjust(x, method="BH", n=length(x))))

##Create the lists for stargazer
fit.list.hav <- list(ha_fit1v, ha_fit2v)
se.list.hav <- list(ha_se1v[,2], ha_se2v[,2])
p.list.hav <- list(ha_pvalues_adjv[,1], ha_pvalues_adjv[,2])
depvar.list.ha <- c("Adult Activity", "Child Activity")
omit.stats <- c("rsq","ser","f")
cov.list <- c("Legal Electrification")

##Create the table
stargazer(title = "Legal Electrification and Household Activity",
          se = se.list.hav,
          header = F,
          fit.list.hav,
          p = p.list.hav,
          covariate.labels = cov.list,
          dep.var.labels = depvar.list.ha,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          model.names = F,
          omit = c("gender","birthplace","age","religion","caste","as.factor"),
          add.lines = list(c("Control Variables", "Yes", "Yes"),
                           c("N Poles", "152", "152"),
                           c("N Villages", "4", "4")),
          notes = "Standard errors clustered at pole level. Pole and village fixed effects included."
)


##TABLE A26##
a_fit1v <- ivreg(appliances ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village), data=analysis)
a_se1v <- coeftest(a_fit1v, cluster.vcov(a_fit1v, as.factor(analysis$Pole)))

#appliance use
a_fit2v <- ivreg(appliance_use ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village), data=analysis)
a_se2v <- coeftest(a_fit2v, cluster.vcov(a_fit2v, as.factor(analysis$Pole)))

##Adjust p values for multiple comparison
a_pvaluesv <- cbind(a_se1v[,4], a_se2v[,4])

#p.adjust(e_pvalues[1,], method="BH", n=length(e_pvalues[1,]))

a_pvalues_adjv <- t(apply(a_pvaluesv, MARGIN=1, function(x) p.adjust(x, method="BH", n=length(x))))

##Create the lists for stargazer
fit.list.av <- list(a_fit1v, a_fit2v)
se.list.av <- list(a_se1v[,2], a_se2v[,2])
p.list.av <- list(a_pvalues_adjv[,1], a_pvalues_adjv[,2])
depvar.list.a <- c("Number of Appliances", "Appliance Use")
omit.stats <- c("rsq","ser","f")
cov.list <- c("Legal Electrification")

##Create the table
stargazer(title = "Legal Electrification and Appliances",
          se = se.list.av,
          header = F,
          fit.list.av,
          p = p.list.av,
          covariate.labels = cov.list,
          dep.var.labels = depvar.list.a,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          model.names = F,
          omit = c("gender","birthplace","age","religion","caste","as.factor"),
          add.lines = list(c("Control Variables", "Yes", "Yes"),
                           c("N Poles", "152", "152"),
                           c("N Villages", "4", "4")),
           
          notes = "Standard errors clustered at pole level. Pole and village fixed effects included."
)


##TABLE A27##
#reliability
sa_fit1v <- ivreg(satisfaction_reliability ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village), data=analysis)
sa_se1v <- coeftest(sa_fit1v, cluster.vcov(sa_fit1v, as.factor(analysis$Pole)))

#cost
sa_fit2v <- ivreg(satisfaction_cost ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village), data=analysis)
sa_se2v <- coeftest(sa_fit2v, cluster.vcov(sa_fit2v, as.factor(analysis$Pole)))

#safety
sa_fit3v <- ivreg(satisfaction_safety ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village), data=analysis)
sa_se3v <- coeftest(sa_fit3v, cluster.vcov(sa_fit3v, as.factor(analysis$Pole)))

#brightness
sa_fit4v <- ivreg(satisfaction_brightness ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village), data=analysis)
sa_se4v <- coeftest(sa_fit4v, cluster.vcov(sa_fit4v, as.factor(analysis$Pole)))


#overall satisfaction
sa_fit5v <- ivreg(satisfaction ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village), data=analysis)
sa_se5v <- coeftest(sa_fit5v, cluster.vcov(sa_fit5v, as.factor(analysis$Pole)))

#change in satisfaction
sa_fit6v <- ivreg(satisfaction_chng ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village), data=analysis)
sa_se6v <- coeftest(sa_fit6v, cluster.vcov(sa_fit6v, as.factor(analysis$Pole)))

#electricity value
sa_fit7v <- ivreg(elec_value ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village), data=analysis)
sa_se7v <- coeftest(sa_fit7v, cluster.vcov(sa_fit7v, as.factor(analysis$Pole)))

#satisfaction with business
sa_fit8v <- ivreg(satisfaction_business ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village), data=analysis)
sa_se8v <- coeftest(sa_fit8v, cluster.vcov(sa_fit8v, as.factor(analysis$Pole)))

#increase in business income
sa_fit9v <- ivreg(income_increase ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village), data=analysis)
sa_se9v <- coeftest(sa_fit9v, cluster.vcov(sa_fit9v, as.factor(analysis$Pole)))

#interest in business
sa_fit10v <- ivreg(business_interest ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village), data=analysis)
sa_se10v <- coeftest(sa_fit10v, cluster.vcov(sa_fit10v, as.factor(analysis$Pole)))

#aspirations/modernity
sa_fit11v <- ivreg(aspirations ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village), data=analysis)
sa_se11v <- coeftest(sa_fit11v, cluster.vcov(sa_fit11v, as.factor(analysis$Pole)))


##Adjust P values
sa_pvaluesv <- cbind(sa_se1v[,4], sa_se2v[,4], sa_se3v[,4], sa_se4v[,4], sa_se5v[,4], sa_se6v[,4], sa_se7v[,4], sa_se8v[,4], sa_se9v[,4], 
                     sa_se10v[,4], sa_se11v[,4])

#p.adjust(e_pvalues[1,], method="BH", n=length(e_pvalues[1,]))

sa_pvalues_adjv <- t(apply(sa_pvaluesv, MARGIN=1, function(x) p.adjust(x, method="BH", n=length(x))))




##Create the lists for stargazer

#Table sa1
fit.list.sa1v <- list(sa_fit1v, sa_fit2v, sa_fit3v, sa_fit4v)
se.list.sa1v <- list(sa_se1v[,2], sa_se2v[,2], sa_se3v[,2], sa_se4v[,2])
p.list.sa1v <- list(sa_pvalues_adjv[,1], sa_pvalues_adjv[,2], sa_pvalues_adjv[,3], sa_pvalues_adjv[,4])
depvar.list.sa1 <- c("Reliability", "Cost", "Safety", "Brightness")

#Create the table
stargazer(title = "Legal Electrification and Satisfaction and Attitudes",
          se = se.list.sa1v,
          header = F,
          fit.list.sa1v,
          p = p.list.sa1v,
          covariate.labels = cov.list,
          dep.var.labels = depvar.list.sa1,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          model.names = F,
          omit = c("gender","birthplace","age","religion","caste","as.factor"),
          add.lines = list(c("Control Variables", "Yes", "Yes", "Yes", "Yes"),
                           c("N Poles", "152", "152", "152", "152"),
                           c("N Villages", "4", "4", "4", "4")),
           
          notes = "Standard errors clustered at pole level. Pole and village fixed effects included."
)

##TABLE A28##
fit.list.sa2v <- list(sa_fit5v, sa_fit6v, sa_fit7v)
se.list.sa2v <- list(sa_se5v[,2], sa_se6v[,2], sa_se7v[,2])
p.list.sa2v <- list(sa_pvalues_adjv[,5], sa_pvalues_adjv[,6], sa_pvalues_adjv[,7])
depvar.list.sa2 <- c("Overall Satisfaction", "Change in Satisfaction", "Electricity Value")

#Create the table
stargazer(title = "Legal Electrification and Satisfaction and Attitudes",
          se = se.list.sa2v,
          header = F,
          fit.list.sa2v,
          p = p.list.sa2v,
          covariate.labels = cov.list,
          dep.var.labels = depvar.list.sa2,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          model.names = F,
          omit = c("gender","birthplace","age","religion","caste","as.factor"),
          add.lines = list(c("Control Variables", "Yes", "Yes", "Yes"),
                           c("N Poles", "152", "152", "152"),
                           c("N Villages", "4", "4", "4")),
           
          notes = "Standard errors clustered at pole level. Pole and village fixed effects included."
)



##TABLE A29##
fit.list.sa3v <- list(sa_fit9v, sa_fit10v, sa_fit11v)
se.list.sa3v <- list(sa_se9v[,2], sa_se10v[,2], sa_se11v[,2])
p.list.sa3v <- list(sa_pvalues_adjv[,9], sa_pvalues_adjv[,10], sa_pvalues_adjv[,11])
depvar.list.sa3 <- c("Increase in Income", "Business Aspiration", "General Aspirations")

#Create the table
stargazer(title = "Legal Electrification and Satisfaction and Attitudes",
          se = se.list.sa3v,
          header = F,
          fit.list.sa3v,
          p = p.list.sa3v,
          covariate.labels = cov.list,
          dep.var.labels = depvar.list.sa3,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          model.names = F,
          omit = c("gender","birthplace","age","religion","caste","as.factor"),
          add.lines = list(c("Control Variables", "Yes", "Yes", "Yes"),
                           c("N Poles", "152", "152", "152"),
                           c("N Villages", "4", "4", "4")),
           
          notes = "Standard errors clustered at pole level. Pole and fixed effects included."
)


##TABLE A30
k_fit1v <- ivreg(knowledge ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole)+as.factor(village), data=analysis)
k_se1v <- coeftest(k_fit1v, cluster.vcov(k_fit1v, as.factor(analysis$Pole)))

#no need for padjust- only 1 model

##Create the lists for stargazer
fit.list.kv <- list(k_fit1v)
se.list.kv <- list(k_se1v[,2])
omit.stats <- c("rsq","ser","f")
cov.list <- c("Legal Electrification")


##Create the table
stargazer(title = "Legal Electrification and Knowledge",
          se = se.list.kv,
          header = F,
          fit.list.kv,
          covariate.labels = cov.list,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          model.names = F,
          omit = c("gender","birthplace","age","religion","caste","as.factor"),
          add.lines = list(c("Control Variables", "Yes"),
                           c("N Poles", "152"),
                           c("N Villages", "4")),
          dep.var.labels = "Knowledge",
          notes = "Standard errors clustered at pole level. Pole and village fixed effects included."
)

##TABLE A31##

#Total Expenditure
e_fit1c <- ivreg(total_expenditure ~ treat + gender + birthplace + age + religion + caste + education + homeBuilt + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + education + homeBuilt +  as.factor(Pole), data=analysis)
e_se1c <- coeftest(e_fit1c, cluster.vcov(e_fit1c, as.factor(analysis$Pole)))

#Food
e_fit2c <- ivreg(food_expenditure ~ treat + gender + birthplace + age + religion + caste + education + homeBuilt + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + education + homeBuilt +  as.factor(Pole), data=analysis)
e_se2c <- coeftest(e_fit2c, cluster.vcov(e_fit2c, as.factor(analysis$Pole)))

#Education
e_fit3c <- ivreg(education_expenditure ~ treat + gender + birthplace + age + religion + caste + education + homeBuilt + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + education + homeBuilt +  as.factor(Pole), data=analysis)
e_se3c <- coeftest(e_fit3c, cluster.vcov(e_fit3c, as.factor(analysis$Pole)))

#kerosene
e_fit4c <- ivreg(kerosene_expenditure ~ treat + gender + birthplace + age + religion + caste + education + homeBuilt + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + education + homeBuilt +  as.factor(Pole), data=analysis)
e_se4c <- coeftest(e_fit4c, cluster.vcov(e_fit4c, as.factor(analysis$Pole)))

#electricity
e_fit5c <- ivreg(electricity_expenditure ~ treat + gender + birthplace + age + religion + caste + education + homeBuilt + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + education + homeBuilt +  as.factor(Pole), data=analysis)
e_se5c <- coeftest(e_fit5c, cluster.vcov(e_fit5c, as.factor(analysis$Pole)))


##Adjust p values for multiple comparison
e_pvaluesc <- cbind(e_se1c[,4], e_se2c[,4], e_se3c[,4], e_se4c[,4])

#p.adjust(e_pvalues[1,], method="BH", n=length(e_pvalues[1,]))

e_pvalues_adjc <- t(apply(e_pvaluesc, MARGIN=1, function(x) p.adjust(x, method="BH", n=length(x))))


##Create the lists for stargazer
fit.list.ec <- list(e_fit1c, e_fit2c, e_fit3c, e_fit4c)
se.list.ec <- list(e_se1c[,2], e_se2c[,2], e_se3c[,2], e_se4c[,2])
p.list.ec <- list(e_pvalues_adjc[,1], e_pvalues_adjc[,2], e_pvalues_adjc[,3], e_pvalues_adjc[,4])
depvar.list.ec <- c("Total Expenditure", "Food Expenditure", "Education Expenditure", "Kerosene Expenditure")
omit.stats <- c("rsq","ser","f")
cov.list <- c("Legal Electrification")

##Create the table
stargazer(title = "Legal Electrification and Expenditure",
          se = se.list.ec,
          header = F,
          fit.list.ec,
          p = p.list.ec,
          covariate.labels = cov.list,
          dep.var.labels = depvar.list.ec,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          model.names = F,
          omit = c("gender","birthplace","age","religion","caste","as.factor", "education", "homeBuilt"),
          add.lines = list(c("Control Variables", "Yes", "Yes", "Yes", "Yes"),
                           c("Control for Home Age and Education", "Yes", "Yes", "Yes", "Yes"),
                           c("N Poles", "152", "152", "152", "152")),
           
          notes = "Standard errors clustered at pole level. Pole fixed effects included."
)

##TABLE A32##
#lamp binary
ke_fit1c <- ivreg(kerosene_lamps ~  treat + gender + birthplace + age + religion + caste + education + homeBuilt + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + education + homeBuilt +  as.factor(Pole), data=analysis)
ke_se1c <- coeftest(ke_fit1c, cluster.vcov(ke_fit1c, as.factor(analysis$Pole)))

#num kerosene lamps
ke_fit2c <- ivreg(num_kerosene_lamps ~  treat + gender + birthplace + age + religion + caste + education + homeBuilt + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + education + homeBuilt +  as.factor(Pole), data=analysis)
ke_se2c <- coeftest(ke_fit2c, cluster.vcov(ke_fit2c, as.factor(analysis$Pole)))

#kerosene lamp hours
ke_fit3c <- ivreg(kerosene_lamp_hours ~  treat + gender + birthplace + age + religion + caste + education + homeBuilt + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + education + homeBuilt +  as.factor(Pole), data=analysis)
ke_se3c <- coeftest(ke_fit3c, cluster.vcov(ke_fit3c, as.factor(analysis$Pole)))

#Other use of Kerosene
ke_fit4c <- ivreg(kerosene_other ~  treat + gender + birthplace + age + religion + caste + education + homeBuilt + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + education + homeBuilt +  as.factor(Pole), data=analysis)
ke_se4c <- coeftest(ke_fit4c, cluster.vcov(ke_fit4c, as.factor(analysis$Pole)))


##Adjust p values for multiple comparison
ke_pvaluesc <- cbind(ke_se1c[,4], ke_se2c[,4], ke_se3c[,4], ke_se4c[,4])

#p.adjust(e_pvalues[1,], method="BH", n=length(e_pvalues[1,]))

ke_pvalues_adjc <- t(apply(ke_pvaluesc, MARGIN=1, function(x) p.adjust(x, method="BH", n=length(x))))


##Create the lists for stargazer
fit.list.kec <- list(ke_fit1c, ke_fit2c, ke_fit3c, ke_fit4c)
se.list.kec <- list(ke_se1c[,2], ke_se2c[,2], ke_se3c[,2], ke_se4c[,2])
p.list.kec <- list(ke_pvalues_adjc[,1], ke_pvalues_adjc[,2], ke_pvalues_adjc[,3], ke_pvalues_adjc[,4])
depvar.list.kec <- c("Kerosene Lamp Use", "Number of Kerosene Lamps", "Kerosene Lamp Hours", "Other Use of Kerosene")
omit.stats <- c("rsq","ser","f")
cov.list <- c("Legal Electrification")

##Create the table
stargazer(title = "Legal Electrification and Kerosene Use",
          se = se.list.kec,
          header = F,
          fit.list.kec,
          p = p.list.kec,
          covariate.labels = cov.list,
          dep.var.labels = depvar.list.kec,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          model.names = F,
          omit = c("gender","birthplace","age","religion","caste","as.factor", "education", "homeBuilt"),
          add.lines = list(c("Control Variables", "Yes", "Yes", "Yes", "Yes"),
                           c("Control for Home Age and Education", "Yes", "Yes", "Yes", "Yes"),
                           c("N Poles", "152", "152", "152", "152")), #add control line
          
          notes = "Standard errors clustered at pole level. Pole fixed effects included."
)

##TABLE A33##
#Lighting hours
l_fit1c <- ivreg(lighting_hours ~  treat + gender + birthplace + age + religion + caste + education + homeBuilt + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + education + homeBuilt +  as.factor(Pole), data=analysis)
l_se1c <- coeftest(l_fit1c, cluster.vcov(l_fit1c, as.factor(analysis$Pole)))

#child lighting hours
l_fit2c <- ivreg(child_lighting ~  treat + gender + birthplace + age + religion + caste + education + homeBuilt + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + education + homeBuilt +  as.factor(Pole), data=analysis)
l_se2c <- coeftest(l_fit2c, cluster.vcov(l_fit2c, as.factor(analysis$Pole)))

#adult lighting hours
l_fit3c <- ivreg(adult_lighting ~  treat + gender + birthplace + age + religion + caste + education + homeBuilt + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + education + homeBuilt +  as.factor(Pole), data=analysis)
l_se3c <- coeftest(l_fit3c, cluster.vcov(l_fit3c, as.factor(analysis$Pole)))


##Adjust p values for multiple comparison
l_pvaluesc <- cbind(l_se1c[,4], l_se2c[,4], l_se3c[,4])

#p.adjust(e_pvalues[1,], method="BH", n=length(e_pvalues[1,]))

l_pvalues_adjc <- t(apply(l_pvaluesc, MARGIN=1, function(x) p.adjust(x, method="BH", n=length(x))))

##Create the lists for stargazer
fit.list.lc <- list(l_fit1c, l_fit2c, l_fit3c)
se.list.lc <- list(l_se1c[,2], l_se2c[,2], l_se3c[,2])
p.list.lc <- list(l_pvalues_adjc[,1], l_pvalues_adjc[,2], l_pvalues_adjc[,3])
depvar.list.lc <- c("Household Light", "Child Light", "Adult Light")
omit.stats <- c("rsq","ser","f")
cov.list <- c("Legal Electrification")

##Create the table
stargazer(title = "Legal Electrification and Lighting Use",
          se = se.list.lc,
          header = F,
          fit.list.lc,
          p = p.list.lc,
          covariate.labels = cov.list,
          dep.var.labels = depvar.list.lc,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          model.names = F,
          omit = c("gender","birthplace","age","religion","caste","as.factor", "education", "homeBuilt"),
          add.lines = list(c("Control Variables", "Yes", "Yes", "Yes"),
                           c("Control for Home Age and Education", "Yes", "Yes", "Yes"),
                           c("N Poles", "152", "152", "152")), #add control list
           
          notes = "Standard errors clustered at pole level. Pole fixed effects included."
)

##TABLE A34##
#adult activity
ha_fit1c <- ivreg(adult_activity ~  treat + gender + birthplace + age + religion + caste + education + homeBuilt + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + education + homeBuilt +  as.factor(Pole), data=analysis)
ha_se1c <- coeftest(ha_fit1c, cluster.vcov(ha_fit1c, as.factor(analysis$Pole)))

#child activity
ha_fit2c <- ivreg(child_activity ~  treat + gender + birthplace + age + religion + caste + education + homeBuilt + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + education + homeBuilt +  as.factor(Pole), data=analysis)
ha_se2c <- coeftest(ha_fit2c, cluster.vcov(ha_fit2c, as.factor(analysis$Pole)))

##Adjust p values for multiple comparison
ha_pvaluesc <- cbind(ha_se1c[,4], ha_se2c[,4])

#p.adjust(e_pvalues[1,], method="BH", n=length(e_pvalues[1,]))

ha_pvalues_adjc <- t(apply(ha_pvaluesc, MARGIN=1, function(x) p.adjust(x, method="BH", n=length(x))))

##Create the lists for stargazer
fit.list.hac <- list(ha_fit1c, ha_fit2c)
se.list.hac <- list(ha_se1c[,2], ha_se2c[,2])
p.list.hac <- list(ha_pvalues_adjc[,1], ha_pvalues_adjc[,2])
depvar.list.hac <- c("Adult Activity", "Child Activity")
omit.stats <- c("rsq","ser","f")
cov.list <- c("Legal Electrification")

##Create the table
stargazer(title = "Legal Electrification and Household Activity",
          se = se.list.hac,
          header = F,
          fit.list.hac,
          p = p.list.hac,
          covariate.labels = cov.list,
          dep.var.labels = depvar.list.hac,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          model.names = F,
          omit = c("gender","birthplace","age","religion","caste","as.factor", "education", "homeBuilt"),
          add.lines = list(c("Control Variables", "Yes", "Yes"),
                           c("Control for Home Age and Education", "Yes", "Yes"),
                           c("N Poles", "152", "152")),
          notes = "Standard errors clustered at pole level. Pole fixed effects included."
)

##TABLE A35##
a_fit1c <- ivreg(appliances ~  treat + gender + birthplace + age + religion + caste + education + homeBuilt + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + education + homeBuilt +  as.factor(Pole), data=analysis)
a_se1c <- coeftest(a_fit1c, cluster.vcov(a_fit1c, as.factor(analysis$Pole)))

#appliance use
a_fit2c <- ivreg(appliance_use ~  treat + gender + birthplace + age + religion + caste + education + homeBuilt + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + education + homeBuilt +  as.factor(Pole), data=analysis)
a_se2c <- coeftest(a_fit2c, cluster.vcov(a_fit2c, as.factor(analysis$Pole)))

##Adjust p values for multiple comparison
a_pvaluesc <- cbind(a_se1c[,4], a_se2c[,4])

#p.adjust(e_pvalues[1,], method="BH", n=length(e_pvalues[1,]))

a_pvalues_adjc <- t(apply(a_pvaluesc, MARGIN=1, function(x) p.adjust(x, method="BH", n=length(x))))

##Create the lists for stargazer
fit.list.ac <- list(a_fit1c, a_fit2c)
se.list.ac <- list(a_se1c[,2], a_se2c[,2])
p.list.ac <- list(a_pvalues_adjc[,1], a_pvalues_adjc[,2])
depvar.list.ac <- c("Number of Appliances", "Appliance Use")
omit.stats <- c("rsq","ser","f")
cov.list <- c("Legal Electrification")

##Create the table
stargazer(title = "Legal Electrification and Appliances",
          se = se.list.ac,
          header = F,
          fit.list.ac,
          p = p.list.ac,
          covariate.labels = cov.list,
          dep.var.labels = depvar.list.ac,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          model.names = F,
          omit = c("gender","birthplace","age","religion","caste","as.factor", "education", "homeBuilt"),
          add.lines = list(c("Control Variables", "Yes", "Yes"),
                           c("Control for Home Age and Education", "Yes", "Yes"),
                           c("N Poles", "152", "152")),
           
          notes = "Standard errors clustered at pole level. Pole fixed effects included."
)

##TABLE A36##
#reliability
sa_fit1c <- ivreg(satisfaction_reliability ~  treat + gender + birthplace + age + religion + caste + education + homeBuilt + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + education + homeBuilt +  as.factor(Pole), data=analysis)
sa_se1c <- coeftest(sa_fit1c, cluster.vcov(sa_fit1c, as.factor(analysis$Pole)))

#cost
sa_fit2c <- ivreg(satisfaction_cost ~  treat + gender + birthplace + age + religion + caste + education + homeBuilt + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + education + homeBuilt +  as.factor(Pole), data=analysis)
sa_se2c <- coeftest(sa_fit2c, cluster.vcov(sa_fit2c, as.factor(analysis$Pole)))

#safety
sa_fit3c <- ivreg(satisfaction_safety ~  treat + gender + birthplace + age + religion + caste + education + homeBuilt + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + education + homeBuilt +  as.factor(Pole), data=analysis)
sa_se3c <- coeftest(sa_fit3c, cluster.vcov(sa_fit3c, as.factor(analysis$Pole)))

#brightness
sa_fit4c <- ivreg(satisfaction_brightness ~  treat + gender + birthplace + age + religion + caste + education + homeBuilt + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + education + homeBuilt +  as.factor(Pole), data=analysis)
sa_se4c <- coeftest(sa_fit4c, cluster.vcov(sa_fit4c, as.factor(analysis$Pole)))


#overall satisfaction
sa_fit5c <- ivreg(satisfaction ~  treat + gender + birthplace + age + religion + caste + education + homeBuilt + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + education + homeBuilt +  as.factor(Pole), data=analysis)
sa_se5c <- coeftest(sa_fit5c, cluster.vcov(sa_fit5c, as.factor(analysis$Pole)))

#change in satisfaction
sa_fit6c <- ivreg(satisfaction_chng ~  treat + gender + birthplace + age + religion + caste + education + homeBuilt + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + education + homeBuilt +  as.factor(Pole), data=analysis)
sa_se6c <- coeftest(sa_fit6c, cluster.vcov(sa_fit6c, as.factor(analysis$Pole)))

#electricity value
sa_fit7c <- ivreg(elec_value ~  treat + gender + birthplace + age + religion + caste + education + homeBuilt + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + education + homeBuilt +  as.factor(Pole), data=analysis)
sa_se7c <- coeftest(sa_fit7c, cluster.vcov(sa_fit7c, as.factor(analysis$Pole)))

#satisfaction with business
sa_fit8c <- ivreg(satisfaction_business ~ treat + gender + birthplace + age + religion + caste + education + homeBuilt + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + education + homeBuilt +  as.factor(Pole), data=analysis)
sa_se8c <- coeftest(sa_fit8c, cluster.vcov(sa_fit8c, as.factor(analysis$Pole)))

#increase in business income
sa_fit9c <- ivreg(income_increase ~  treat + gender + birthplace + age + religion + caste + education + homeBuilt + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + education + homeBuilt +  as.factor(Pole), data=analysis)
sa_se9c <- coeftest(sa_fit9c, cluster.vcov(sa_fit9c, as.factor(analysis$Pole)))

#interest in business
sa_fit10c <- ivreg(business_interest ~  treat + gender + birthplace + age + religion + caste + education + homeBuilt + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + education + homeBuilt +  as.factor(Pole), data=analysis)
sa_se10c <- coeftest(sa_fit10c, cluster.vcov(sa_fit10c, as.factor(analysis$Pole)))

#aspirations/modernity
sa_fit11c <- ivreg(aspirations ~  treat + gender + birthplace + age + religion + caste + education + homeBuilt + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + education + homeBuilt +  as.factor(Pole), data=analysis)
sa_se11c <- coeftest(sa_fit11c, cluster.vcov(sa_fit11c, as.factor(analysis$Pole)))


##Adjust P values
sa_pvaluesc <- cbind(sa_se1c[,4], sa_se2c[,4], sa_se3c[,4], sa_se4c[,4], sa_se5c[,4], sa_se6c[,4], sa_se7c[,4], sa_se8c[,4], sa_se9c[,4], 
                     sa_se10c[,4], sa_se11c[,4])

#p.adjust(e_pvalues[1,], method="BH", n=length(e_pvalues[1,]))

sa_pvalues_adjc <- t(apply(sa_pvaluesc, MARGIN=1, function(x) p.adjust(x, method="BH", n=length(x))))




##Create the lists for stargazer


fit.list.sa1c <- list(sa_fit1c, sa_fit2c, sa_fit3c, sa_fit4c)
se.list.sa1c <- list(sa_se1c[,2], sa_se2c[,2], sa_se3c[,2], sa_se4c[,2])
p.list.sa1c <- list(sa_pvalues_adjc[,1], sa_pvalues_adjc[,2], sa_pvalues_adjc[,3], sa_pvalues_adjc[,4])
depvar.list.sa1c <- c("Reliability", "Cost", "Safety", "Brightness")
cov.list <- c("Legal Electrification")

#Create the table
stargazer(title = "Legal Electrification and Satisfaction and Attitudes",
          se = se.list.sa1c,
          header = F,
          fit.list.sa1c,
          p = p.list.sa1c,
          covariate.labels = cov.list,
          dep.var.labels = depvar.list.sa1c,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          model.names = F,
          omit = c("gender","birthplace","age","religion","caste","as.factor", "education", "homeBuilt"),
          add.lines = list(c("Control Variables", "Yes", "Yes", "Yes", "Yes"),
                           c("Control for Home Age and Education", "Yes", "Yes", "Yes", "Yes"),
                           c("N Poles", "152", "152", "152", "152")),
           
          notes = "Standard errors clustered at pole level. Pole fixed effects included."
)


##TABLE A37##
fit.list.sa2c <- list(sa_fit5c, sa_fit6c, sa_fit7c)
se.list.sa2c <- list(sa_se5c[,2], sa_se6c[,2], sa_se7c[,2])
p.list.sa2c <- list(sa_pvalues_adjc[,5], sa_pvalues_adjc[,6], sa_pvalues_adjc[,7])
depvar.list.sa2c <- c("Overall Satisfaction", "Change in Satisfaction", "Electricity Value")

#Create the table
stargazer(title = "Legal Electrification and Satisfaction and Attitudes",
          se = se.list.sa2c,
          header = F,
          fit.list.sa2c,
          p = p.list.sa2c,
          covariate.labels = cov.list,
          dep.var.labels = depvar.list.sa2c,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          model.names = F,
          omit = c("gender","birthplace","age","religion","caste","as.factor", "education", "homeBuilt"),
          add.lines = list(c("Control Variables", "Yes", "Yes", "Yes"),
                           c("Control for Home Age and Education", "Yes", "Yes", "Yes"),
                           c("N Poles", "152", "152", "152")),
           
          notes = "Standard errors clustered at pole level. Pole fixed effects included."
)



##TABLE A38##
fit.list.sa3c <- list(sa_fit9c, sa_fit10c, sa_fit11c)
se.list.sa3c <- list(sa_se9c[,2], sa_se10c[,2], sa_se11c[,2])
p.list.sa3c <- list(sa_pvalues_adjc[,9], sa_pvalues_adjc[,10], sa_pvalues_adjc[,11])
depvar.list.sa3c <- c("Increase in Income", "Business Aspiration", "General Aspirations")

#Create the table
stargazer(title = "Legal Electrification and Satisfaction and Attitudes",
          se = se.list.sa3c,
          header = F,
          fit.list.sa3c,
          p = p.list.sa3c,
          covariate.labels = cov.list,
          dep.var.labels = depvar.list.sa3c,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          model.names = F,
          omit = c("gender","birthplace","age","religion","caste","as.factor", "education", "homeBuilt"),
          add.lines = list(c("Control Variables", "Yes", "Yes", "Yes"),
                           c("Control for Home Age and Education", "Yes", "Yes", "Yes"),
                           c("N Poles", "152", "152", "152")),
           
          notes = "Standard errors clustered at pole level. Pole fixed effects included."
)

##TABLE A39##
k_fit1c <- ivreg(knowledge ~  treat + gender + birthplace + age + religion + caste + education + homeBuilt + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + education + homeBuilt +  as.factor(Pole), data=analysis)
k_se1c <- coeftest(k_fit1c, cluster.vcov(k_fit1c, as.factor(analysis$Pole)))

#no need for padjust- only 1 model

##Create the lists for stargazer
fit.list.kc <- list(k_fit1c)
se.list.kc <- list(k_se1c[,2])
omit.stats <- c("rsq","ser","f")
cov.list <- c("Legal Electrification")


##Create the table
stargazer(title = "Legal Electrification and Knowledge",
          se = se.list.kc,
          header = F,
          fit.list.kc,
          covariate.labels = cov.list,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          model.names = F,
          omit = c("gender","birthplace","age","religion","caste","as.factor", "education", "homeBuilt"),
          add.lines = list(c("Control Variables", "Yes"),
                           c("Control for Home Age and Education", "Yes"),
                           c("N Poles", "152")),
          dep.var.labels = "Knowledge",
          notes = "Standard errors clustered at pole level. Pole fixed effects included."
)


##TABLE A40##
#Total Expenditure
e_fit1_log <- ivreg(te_log ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
e_se1_log <- coeftest(e_fit1_log, cluster.vcov(e_fit1_log, as.factor(analysis$Pole)))

#Food
e_fit2_log <- ivreg(fe_log ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
e_se2_log <- coeftest(e_fit2_log, cluster.vcov(e_fit2_log, as.factor(analysis$Pole)))


e_log_pvalues <- cbind(e_se1_log[,4], e_se2_log[,4])


e_log_pvalues_adj <- t(apply(e_log_pvalues, MARGIN=1, function(x) p.adjust(x, method="BH", n=length(x))))



#then redo them with outliers dropped- dropped all above 15000
outlier_te <- analysis[!analysis$total_expenditure>15000,]


outlier_fe <- analysis[!analysis$food_expenditure>15000,]
#Total Expenditure
e_fit1_o <- ivreg(total_expenditure ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=outlier_te)
e_se1_o <- coeftest(e_fit1_o, cluster.vcov(e_fit1_o, as.factor(outlier_te$Pole)))

#Food
e_fit2_o <- ivreg(food_expenditure ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=outlier_fe)
e_se2_o <- coeftest(e_fit2_o, cluster.vcov(e_fit2_o, as.factor(outlier_fe$Pole)))

e_o_pvalues <- cbind(e_se1_o[,4], e_se2_o[,4])


e_o_pvalues_adj <- t(apply(e_o_pvalues, MARGIN=1, function(x) p.adjust(x, method="BH", n=length(x))))


###CREATE TABLE
##Create the lists for stargazer
fit.list.eadj <- list(e_fit1_log, e_fit2_log, e_fit1_o, e_fit2_o)
se.list.eadj <- list(e_se1_log[,2], e_se2_log[,2], e_se1_o[,2], e_se2_o[,2])
p.list.eadj <- list(e_log_pvalues_adj[,1], e_log_pvalues_adj[,2], e_o_pvalues_adj[,1], e_o_pvalues_adj[,2])
depvar.list.eadj <- c("Total Expenditure (log)", "Food Expenditure (log)", "Total Expenditure", "Food Expenditure")
omit.stats <- c("rsq","ser","f")
cov.list <- c("Legal Electrification")

stargazer(title = "Adjusted Legal Electrification and Expenditure",
          se = se.list.eadj,
          header = F,
          fit.list.eadj,
          p = p.list.eadj,
          covariate.labels = cov.list,
          dep.var.labels = depvar.list.eadj,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          model.names = F,
          omit = c("gender","birthplace","age","religion","caste","as.factor"),
          add.lines = list(c("Control Variables", "Yes", "Yes", "Yes", "Yes"),
                           c("N Poles", "152", "152", "152", "152")),
          notes = "Standard errors clustered at pole level. Pole fixed effects included."
)

##TABLE A41##
##Objects come from code under Table A3

fit.list.zsc <- list(ex_zsc_reg, k_zsc_reg, l_zsc_reg, a_zsc_reg, ha_zsc_reg, sa_zsc_reg, k_fit1)
se.list.zsc <- list(ex_zsc_reg_se[,2], k_zsc_reg_se[,2], l_zsc_reg_se[,2], a_zsc_reg_se[,2], ha_zsc_reg_se[,2], sa_zsc_reg_se[,2], k_se1[,2])
p.list.zsc <- list(ex_zsc_reg_se[,4], k_zsc_reg_se[,4], l_zsc_reg_se[,4], a_zsc_reg_se[,4], ha_zsc_reg_se[,4], sa_zsc_reg_se[,4], k_se1[,4])
depvar.list.zsc <- c("Total Expenditure", "Kerosene", "Artificial Lighting", "Appliances", "Household Activity", "Satisfaction and Attitudes", "Knowledge")
omit.stats <- c("rsq","ser","f")
cov.list <- c("Legal Electrification")



##Create the table
stargazer(title = "Legal Electrification and Mean Effects Indices",
          se = se.list.zsc,
          header = F,
          fit.list.zsc,
          p = p.list.zsc,
          covariate.labels = cov.list,
          dep.var.labels = depvar.list.zsc,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          model.names = F,
          omit = c("gender","birthplace","age","religion","caste","as.factor"),
          add.lines = list(c("Control Variables", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("N Poles", "152", "152", "152", "152", "152", "152", "152"),
                           c("N Villages", "4", "4", "4", "4", "4", "4", "4")),
          notes = "Standard errors clustered at pole level. Pole and village fixed effects included."
)


##TABLE A42##
lx1 <- ivreg(lighting_business ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
lxse1 <- coeftest(lx1, cluster.vcov(lx1, as.factor(analysis$Pole)))

lx2 <- ivreg(lighting_friends_family ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
lxse2 <- coeftest(lx2, cluster.vcov(lx2, as.factor(analysis$Pole)))

lx3 <- ivreg(lighting_studying ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
lxse3 <- coeftest(lx3, cluster.vcov(lx3, as.factor(analysis$Pole)))

lx4 <- ivreg(lighting_reading ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
lxse4 <- coeftest(lx4, cluster.vcov(lx4, as.factor(analysis$Pole)))

lx5 <- ivreg(lighting_media ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
lxse5 <- coeftest(lx5, cluster.vcov(lx5, as.factor(analysis$Pole)))

lx6 <- ivreg(lighting_other ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
lxse6 <- coeftest(lx6, cluster.vcov(lx6, as.factor(analysis$Pole)))

##Create the lists for stargazer
fit.list.lx <- list(lx1, lx2, lx3, lx4, lx5, lx6)
se.list.lx <- list(lxse1[,2], lxse2[,2], lxse3[,2], lxse4[,2], lxse5[,2], lxse6[,2])
depvar.list.lx <- c("Business", "Friends/Family", "Studying", "Reading", "Media", "Other")
omit.stats <- c("rsq","ser","f")

##Create the table
stargazer(title = "Decomposed Measures of Artificial Lighting",
          se = se.list.lx,
          header = F,
          fit.list.lx,
          
          covariate.labels = cov.list,
          dep.var.labels = depvar.list.lx,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          model.names = F,
          omit = c("gender","birthplace","age","religion","caste","as.factor"),
          add.lines = list(c("Control Variables", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("N Poles", "152", "152",  "152",  "152",  "152",  "152")),
          notes = "Standard errors clustered at pole level. Pole fixed effects included."
)


##TABLE A43##

ap1 <- ivreg(mobile ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
apse1 <- coeftest(ap1, cluster.vcov(ap1, as.factor(analysis$Pole)))


ap2 <- ivreg(ic_bulb ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
apse2 <- coeftest(ap2, cluster.vcov(ap2, as.factor(analysis$Pole)))


ap3 <- ivreg(cfl_bulb ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
apse3 <- coeftest(ap3, cluster.vcov(ap3, as.factor(analysis$Pole)))



ap4 <- ivreg(led_light ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
apse4 <- coeftest(ap4, cluster.vcov(ap4, as.factor(analysis$Pole)))


ap5 <- ivreg(tube_light ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
apse5 <- coeftest(ap5, cluster.vcov(ap5, as.factor(analysis$Pole)))


ap6 <- ivreg(fans ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
apse6 <- coeftest(ap6, cluster.vcov(ap6, as.factor(analysis$Pole)))


ap7 <- ivreg(iron ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
apse7 <- coeftest(ap7, cluster.vcov(ap7, as.factor(analysis$Pole)))


ap8 <- ivreg(refrigerator ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
apse8 <- coeftest(ap8, cluster.vcov(ap8, as.factor(analysis$Pole)))


ap9 <- ivreg(television ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
apse9 <- coeftest(ap9, cluster.vcov(ap9, as.factor(analysis$Pole)))


ap10 <- ivreg(radio ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
apse10 <- coeftest(ap10, cluster.vcov(ap10, as.factor(analysis$Pole)))


ap11 <- ivreg(cooler ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
apse11 <- coeftest(ap11, cluster.vcov(ap11, as.factor(analysis$Pole)))



ap13 <- ivreg(electric_stove ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
apse13 <- coeftest(ap13, cluster.vcov(ap13, as.factor(analysis$Pole)))



ap15 <- ivreg(water_pump ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
apse15 <- coeftest(ap15, cluster.vcov(ap15, as.factor(analysis$Pole)))


ap16 <- ivreg(other ~ treat + gender + birthplace + age + religion + caste + as.factor(Pole) | forcing + gender + birthplace + age + religion + caste + as.factor(Pole), data=analysis)
apse16 <- coeftest(ap16, cluster.vcov(ap16, as.factor(analysis$Pole)))

fit.list.ap1 <- list(ap1, ap2, ap3, ap4, ap5)
se.list.ap1 <- list(apse1[,2], apse2[,2], apse3[,2], apse4[,2], apse5[,2])
depvar.list.ap1 <- c("Mobile Phones", "Incandescent bulbs", "CFL bulbs", "LED Lights", "Tube Lights")
omit.stats <- c("rsq","ser","f")

##Create the table
stargazer(title = "Decomposed Measures of Appliances",
          se = se.list.ap1,
          header = F,
          fit.list.ap1,
          
          covariate.labels = cov.list,
          dep.var.labels = depvar.list.ap1,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          model.names = F,
          omit = c("gender","birthplace","age","religion","caste","as.factor"),
          add.lines = list(c("Control Variables", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("N Poles", "152", "152",  "152",  "152",  "152",  "152")),
          notes = "Standard errors clustered at pole level. Pole fixed effects included."
)


##TABLE A44##
fit.list.ap2 <- list(ap6, ap7, ap8, ap9, ap10,)
se.list.ap2 <- list(apse6[,2], apse7[,2], apse8[,2], apse9[,2], apse10[,2])
depvar.list.ap2 <- c("Fans", "Electric Irons", "Refrigerators", "Television", "Radios")
omit.stats <- c("rsq","ser","f")

##Create the table
stargazer(title = "Decomposed Measures of Appliances",
          se = se.list.ap2,
          header = F,
          fit.list.ap2,
          
          covariate.labels = cov.list,
          dep.var.labels = depvar.list.ap2,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          model.names = F,
          omit = c("gender","birthplace","age","religion","caste","as.factor"),
          add.lines = list(c("Control Variables", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("N Poles", "152", "152",  "152",  "152",  "152",  "152")),
          notes = "Standard errors clustered at pole level. Pole fixed effects included."
)

##TABLE A45
fit.list.ap3 <- list(ap11, ap13, ap15, ap16)
se.list.ap3 <- list(apse11[,2], apse13[,2], apse15[,2], apse16[,2])
depvar.list.ap3 <- c("Coolers", "Electric Stove", "Electric Water Pumps", "Others")
omit.stats <- c("rsq","ser","f")

##Create the table
stargazer(title = "Decomposed Measures of Appliances",
          se = se.list.ap3,
          header = F,
          fit.list.ap3,
          
          covariate.labels = cov.list,
          dep.var.labels = depvar.list.ap3,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          model.names = F,
          omit = c("gender","birthplace","age","religion","caste","as.factor"),
          add.lines = list(c("Control Variables", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("N Poles", "152", "152",  "152",  "152",  "152",  "152")),
          notes = "Standard errors clustered at pole level. Pole fixed effects included."
)


##TABLE A46: SEE .do##

##TABLE A47: SEE OOS.R##